Este documento es una versión actualizada -en varios sentidos- de unos materiales de análisis de microparrays que escribimos allá por 2010, usando LateX y Sweave, junto con la profesora M. Carme Ruíz de Villa para la asignatura “Análisis de Microarrays” del posgrado de Bioinformática de la Universitat Oberta de Catalunya .
En estos años han cambiado muchas cosas. El posgrado es ahora un máster interuniversitario entre la UOC y la Universidad de Barcelona (UB). Muchos estadísticos i/o bioinformáticos -entre los que me incluyo- han cambiado de Latex/Sweave a RStudio/RMarkdown y Github. Mi compañera de escritura, Mamen, se ha jubilado, pero un nuevo compañero, Ricardo, se ha apuntado a la aventura de revivir estos materiales. Y aunque sigue sin estar claro si los microarrays también se jubilaran pronto, el campo de las ómicas ha crecido de forma explosiva. El “Análisis de datos ómicos”, título de este nuevo documento ya no se centra en los microarrays. También trata de RNA-seq, Proteómica, Metagenómica, Metabolómica o Epigenómica por citar sólo algunos de los más populares.
Es obvio que un curso de “Análisis de datos ómicos” no puede abarcar todas estas disciplinas, aunque también lo es que comparten muchos elementos comunes.
El plan de desarrollo para esta nueva edición, que podra seguirse en el repositorio de Github (), tendrá dos fases:
Previsiblemente este documento, con sus errores e imperfecciones se quedará como “open source” para uso de quien lo desee sin tener que pasar por el doloroso y anquilosante proceso de convertirlo en un libro, que, en este caso, empezaría a quedar obsoleto al minuto de estar acabado.
Este capítulo es una corta transición entre la primera parte del curso, en la que se han presentado los conceptos y herramientas básicos y la segunda en donde se presentan por separado y con mayor detalle los métodos de análisis de datos de microarrays.
Su objetivo por tanto es ofrecer una visión de conjunto que sirva de guía (“roadmap”) para los capítulos siguientes de forma que sin perder el detalle de cada uno de ellos tengamos conciencia de en que punto del proceso general nos encontramos.
El capítulo se estructura en dos partes. En la primera se presentan brevemente algunos de los problemas que típicamente se puede querer estudiar con microarrays u otras técnicas similares de análisis de datos de alto rendimiento. A continuación se presenta lo que se ha llamdo aquí el proceso de análisis de microarrays. Finalmente se intoducen algunos casos reales que, a modo de ejemplo se utilizarán en los capítulos siguientes.
Los microarrays y otras tecnologías de alto rendimiento se han aplicado a multitud de investigaciones, de tipus muy diversos que van desde estudio del cancer (Alizadeh et al. (2000), Golub et al. (1999), van ’t Veer et al. (2002)) al de la germinación y la maduración del tomate (Moore et al. (2002)). A pesar de ello no resulta complicado clasificar los estudios realizados en algunos de los grandes bloques que se describen a continuación. La clasificación está basada en el excelente texto de Simon y colegas (Simon et al. (2003)) y aunque se origina en problemas de microarrays se puede aplicar fácilmente a estudios de genómica o ultrasecuenciación.
El objetivo de los estudios comparativos es determinar si los perfiles de expresión génica difieren entre grupos previamente identificados. También se conoce estos estudios como de selección de genes diferencialmente expresados y son, sin duda los más habituales. Los grupos pueden representar una gran variedad de condiciones, desde distintos tejidos a distintos tratamientos o múltiples combinaciones de factores experimentales.
El análisis de este tipo de experimentos, que se describe en el capítulo sobre selección de genes diferencialmente expresados utiliza herramientas estadísticas como las pruebas de comparación de grupos paramétricas (t de Student) o no (test de Mann-Whitney) o diversos métodos de análisis de la varianza.
Entre los ejemplos de la sección @ref{c4examples} los casos @ref{celltypes}, @ref{estrogen} o @ref{CCL4} hacen referencia a estudios comparativos.
La predicción de clase puede confundirse con la selección de genes en tanto que disponemos de clases predefinidas pero su objetivo es distinto, ya que no pretende simplemente buscar genes cuya expresión sea distinta entre dichos grupos sino genes que puedan ser utilizados para identificar a que clase pertenece un “nuevo” individuo dado cuya clase es “a priori” desconocida. El proceso de predicción suele empezar con una selección de genes informativos, que pueden ser, o no, los mismos que se obtendrían si aplicáramos los métodos del apartado anterior, seguida de la construcción de un modelo de predicción y, lo que es más importante, de la verificación o validación de dicho modelo con unos datos nuevos independientes de los utilizados para el desarrollo del modelo.
Aunque el interés de la predicción de clase es muy alto se trata de un procedimiento mucho más complejo y con más posibilidades de error que la simple selección de genes diferencialmente expresados.
Entre los ejemplos de la sección @ref{c4examples} el caso @ref{golub} trata de un problema de predicción, a la vez que uno de descubrimiento de clases.
Un problema distinto a los descritos se presenta cuando no se conoce las clases en que se agrupan los individuos. En este caso de lo que se trata es de encontrar grupos entre los datos que permitan reunir a los individuos más parecidos entre si y distintos de los de los demás grupos. Los métodos estadísticos que se emplearan en estos casos se conocen como análisis de clusters y no son tan complejos como los de predicción de clase aunque algunos aspectos como por ejemplo la definición del número de grupos no resulta tampoco sencillo.
Entre los ejemplos de la sección tanto el caso golub, en parte, como el @ref{breastTum} tratan problemas de descubrimiento de clases.
Una curiosidad del campo de la estadística es que el término clasificación aparece usado de forma indistinta para referirse a problemas de predicicón de clase o de descubrimiento de clase.
Una vez identificados los principales tipos de estudios quedan muchos que no coinciden plenamente con ninguno de ellos. Sin entrar en detalles podemos señalar los estudios de evolución a lo largo del tiempo (“time course”), los de significación biológica (“Gene Enrichment Analysis”, “Gene Set Enrichment Analysis”, …) , los que buscan relaciones entre los genes (“network analysis” o “pathway analysis”). De momento con conocer e identificar los tres grandes bloques mencionados resultará más que suficiente.
Una de las dificultades con que se encuentra la persona que comienza en el análisis de datos de microarrays es de donde obtener ejemplos concretos con los que prácticar las técnicas que está aprendiendo.
No es difícil encontrar datos de microarrays en internet por lo que se han seleccionado algunos conjuntos de datos interesantes para utilizarlos de ejemplo a lo largo del curso. Algunos de éstos son “populares” en el sentido de que han sido utilizados en diversas ocasiones y por lo tanto se encuentran bien documentados. Otros se han escogido simplemente porque ilustran bien algunas de las ideas que se desea exponer o por su accesibilidad.
Todos los datos corresponden a investigaciones publicadas por lo que no se describen exhaustivamente sinó que se expone brevemente el origen y objetivos del trabajo –incluyendo su clasificación según los grupos definidos en la sección anterior– y las características perinentes para el análisis como el tipo de microarrays, los grupos –si los hay– o como acceder a los datos.
Este conjunto de datos, que se denominará “celltypes”, corresponde a un estudio realizado por Chelvarajan y sus colegas (Chelvarajan et al. (2006)) que analizaron el efecto de la estimulación con lipopolisacáridos en la regulación por parte de citoquinas de ciertos procesos biológicos relacionados con la inflamación.
Este estudio es del tipo “class comparison” es decir su principal objetivo es la obtención de genes diferencialmente expresados entre dos o más condiciones.
Los datos se encuentran disponibles en la base de datos pública GEne Expression Omnibus mantenida por el National Institute of Health (NIH), pero pueden descargarse de la página de materiales del curso para garantizar su disponibilidad.
A finales de los años 90, Todd Golub y sus colaboradores (Golub et al. (1999)) realizaron uno de los estudios más populares hasta el momento con datos de microarrays. En él utilizaron microarrays de oligonucleótidos para 6817 genes humanos para mirar de encontrar una forma de distinguir (clasificar) tumores de pacientes con leucemia linfoblástica aguda (ALL) de aquellos que sufrían de leucemia mieloide aguda (AML). Además se interesaba por la posibilidad de descubrir subgrupos de forma que pudieran definirse variantes de cada una de estas patologías a nivel molecular.
La diversidad de objetivos del estudio lleva a clasificarlo tanto entre los del tipo de predicción de clase como entre los que buscan descubrir nuievas clases o grupos en los datos.
Los datos de este estudio se encuentran disponibles en la web del instituto Broad, en donde se llevó a cabo (http://www.broadinstitute.org). También se encuentra disponible un paquete de R denominado ALL que permite utilizarlos directamente usando R y Bioconductor.
Scholtens y colegas (Scholtens et al. (2004)) describen un estudio sobre el efecto de un tratamiento con estrógenos y del tiempo transcurrido desde el tratamiento en la expresión génica en pacientes de cáncer de mama. Los investigadores supusieron que los genes asociados con una respuesta temprana podrían considerarse dianas directas del tratamiento, mientras que los que tardaron más en hacerlo podrían considerarse objetivos secundarios correspondientes a dianas más alejadas en las vías metabólicas.
Estos datos han sido utilizados multitud de veces en los cursos de análisis de microarrays realizados por el proyecto Bioconductor y se encuentran disponibles en forma de paquete de R, el paquete estrogen. Una cracaterística importante de este paquetes es el hecho de que en vez de los datos procesados proporciona los datos “crudos”en forma de archivos .CEL de Affymetrix. Esto permite una mayor flexibilidad a la hora de reutilizarlos lo que explica su popularidad.
Holger y colegas de la empresa LGC Ltd. en Teddington, Inglaterra realizaron unos expeimentos con microarrays de dos colores en los que trataron hepatocitos de ratón con tetraclorido de carbono (CCL4) o con dimetilsulfóxido (DMSO). El tetraclorido de carbono fue ampliamente utilizado en productos de limpieza o refrigeración para el hogar hasta que se detectó que podía tener efectos tóxicos e incluso cancerígenos. El DMSO es un solvente similar, sin efectos tóxicos conocidos, que se utilizó como control negativo.
Los datos de este estudio no han sido publicados pero se encuentran disponibles en el paquete CCL4 de bioconductor. Su interés reside por un lado en que se trata de datos de microarrays de dos colores de la marca Agilent –en un momento en que la mayoría de estudios se realizan con datos de un color. Aparte de esto cabe resaltar el hecho de que el paquete incluye, de forma similar al anterior, los datos “crudos” en forma de archivos de tipo “Genepix” uno de los programas populares para escanear imágenes generadas con microarrays de dos colores.
Este estudio es también un estudio de comparación de clases cuyo objetivo principal es la selección de genes cuya expresión se asocia al tratamiento con CCL4 o DMSO.
Los datos de este ejemplo denominado kidney son datos ya normalizados referidos a la expresión de los genes en distintos momentos del ciclo delular de la levadura e decir desde que concluye la división celular hasta que se inicia la siguiente.
Los datos puede descargarse desde la página del proyecto “Yeast Cell Cycle Project” (Proyecto de estudio del ciclo celular de la levadura) en la dirección:
La tabla resume la lista de ejemplos que se utilizan en este manual indicando el nombre con que nos referiremos de aquí en adelante a cada conjunto de datos as' como algunas de sus carcaterísticas.
La tabla siguiente resume los “datasets” utilizados en este manual. Aparte del nombre (arbitrario y “mnemotécnico”) se indica el tipo de microarrays, el número de muestras, y el tipo de problema para el que se utilizaron originalmente.
| Nombre | Tipo (Modelo) | N. Muestras | Tipo de estudio |
|---|---|---|---|
celltypes |
Un color (Affy, Mouse 4302) | 12 | Comparativo |
golub |
Un color (Affy, HGU95A) | 38 | Clasificación |
estrogen |
Un color (Affy, HGU95A) | 8 | Comparativo |
CCL4 |
Dos colores (Agilent, WG Rat Microarray) | 16 | Comparativo |
breastTum |
Un color (Affy, HGU95A) | 49 | Clasificación |
Una vez descritos los tipos de análisis y algunos ejemplos podemos pasar a describir el proceso de análisis de microarrays que se resume brevemente en la figura @ref{c04analysisProcess}.
El análisis de microarrays, como la mayoría de análisis debe proceder de forma ordenada y siguiendo el método científico:
Tal como ilustra la figura @ref(fig:c04analysisProcess>) el análisis de microarrays puede ser fácilmente visualizado como un proceso que empieza por una pregunta biológica y concluye con una interpretación de los resultados de los análisis que, de alguna forma, confiamos nos acerque un poco a la respuesta de la pregunta inicial.
Figure 2.1: Proceso de análisis de microarrays
El proceso descrito es básicamente una forma razonable de proceder en general. Los microarrays y otros datos genómicos son diferentes en su naturaleza de los datos clásicos alrededor de los que se han desarrollado la mayor parte de técnicas estadísticas. En consecuencia, en muchos casos ha sido necesario adaptar las técnicas existentes o desarrollar otras nuevas para adecuarse a las nuevas situaciones encontradas. Esto ha determinado que existan muchos métodos para cada una de los pasos descritos anteriormente lo que da lugar a una grandísima cantidad de posibilidades.
En la práctica lo que suele hacerse es optar por utilizar algunos de los métodos en los que hay un cierto consenso acerca de su calidad y utilidad para cada problema. Allison (Allison et al. (2006)) repasa los puntos principales de este consenso dando una lista de puntos a tener en cuenta en cualquier estudio que utilice microarrays. Imbeaud y Auffray (Imbeaud and Auffray (2005)) citan una lista de hasta 39 puntos que uno debe seguir en un experimento con microarrays para usar “buenas prácticas”.
Finalmente Zhu y otros (Zhu, Miecznikowski, and Halfon (2010)) utilizan un conjunto de arrays con valores de expresión conocidos para proponer los que, a su parecer, resultan los métodos más apropiados para cada etapa desde la corrección del backround hasta la selección de genes diferencialmente expresados. La figura 2.2 ilustra algunas de las opciones sugeridas por dichos autores.
Figure 2.2: Diseño de arrays
Los capítulos que siguen al presente proceden aproximadamente en el orden del proceso descrito en 2.1.
En este capítulo se examinan unas componentes clave del análisis de microarrays el diseño de experimentos que, no tan solo es crucial para una buena recogida de información, sino que marca todo el proceso desde el preprocesado al análisis final.
Los datos genómicos son muy variables. La figura 3.1, tomada de (???) ilustra algunas posibles fuentes de variabilidad, desde que se inicia el experimento hasta que se lee la información con el escáner. Sin tener que entrar en cómo influye cada una de estas fuents o cuan grande es su influencia lo que muestra esta figura es que este tipo de experimentos se enfrenta a múltiples fuentes de error, por lo que puede beneficiarse de un correcto diseño de experimentos.
Figure 3.1: Proceso de análisis de microarrays
Habitualmente, en la mayor parte de situaciones experimentales, podemos distinguir entre variabilidad sistemática y variabilidad aleatoria.
La variabilidad sistemática es principalmente debida a procesos técnicos mientras que la variabilidad aleatoria es atribuible tanto a razones técnicas como biológicas. Podemos encontrar ejemplos de variabilidad sistemática en la extracción del ARN, marcaje o foto detección. La variabilidad aleatoria puede estar relacionada con muchos factores tales como la calidad del ADN o las características biológicas de las muestras.
La manera natural de manejar la variabilidad aleatoria es, por supuesto, la utilización de un diseño experimental apropiado y el apoyo para obtener conclusiones de unas herramientas estadísticas adecuadas. Los problemas relacionados con el c de los experimentos los discutiremos en esta sección y los relacionados con la aplicación de métodos estadísticos serán tratados en la sección métodos estadísticos.
Tradicionalmente, se estiman las correcciones de la variabilidad sistemática a partir de los datos, en lo que se llama genéricamente “calibrado”. En este contexto, hablaremos de “normalización”, que se tratará más adelante, en el capítulo dedicado al preprocesado de los datos.
Empezaremos por definir conceptos que aparecen de forma usual cuando se plantea un diseño:
Al planificar un experimento hay tres principios básicos que se deben tener siempre en cuenta: La replicación, el control local o “bloqueo” y la aleatorización.
Aplicados correctamente dichos principios garantizan diseños experimentales eficientes.
La replicación o repetición de un experimento de forma idéntica en un número determinado de unidades, es la que permite la realización de un posterior análisis estadístico.
La utilización de replicas es importante para incrementar la precisión, obtener suficiente potencia en los tests y como base para los procedimientos de inferencia. Normalmente, distinguimos dos tipos de replicas en el análisis de microarrays:
La replicación técnica se utiliza cuando estamos tratando réplicas del mismo material biológico. Puede ser tanto la replicación de spots en el mismo chip como la de diferentes alícuotas de la misma muestra hibridadas en diferentes microarrays.
La replicación biológica se da cuando se toman medidas sobre muestras independientes, normalmente sobre individuos diferentes.
La replicación técnica permite la estimación del error a nivel de medida, mientras que la replicación biológica permite estimar la variabilidad a nivel de población.
Figure 3.2: Replicas biológicas vs réplicas técnicas
Sorprendentemente, los primeros experimentos de microarrays utilizaban pocas o ninguna replica biológica. La principal explicación para este hecho, además de la falta de conocimiento estadístico, fue el alto coste de cada microarray.
En pocos años, la necesidad de las réplicas ha llegado a ser indiscutible, y al mismo tiempo, el coste de los chips ha disminuido considerablemente. Actualmente, es normal la utilización de, al menos, tres a cinco replicas por condición experimental, aunque a este consenso se ha llegado más por razones empíricas que por la disponibilidad de modelos apropiados para el análisis de la potencia y tamaño de la muestra.
Recientemente, se ha producido una importante afluencia de artículos describiendo métodos para el análisis de la potencia y tamaño de la muestra. A pesar de su variedad, ningún método aparece como candidato claro para ser utilizado en situaciones prácticas. Esto se debe, probablemente, a la complejidad de los datos de microarrays, básicamente por que los genes no son independientes. Por tanto, estas estructuras de correlación existen en los datos, pero la mayor parte de las dependencias son desconocidas por lo que la estimación de estas estructuras es muy complicada.
Como indicaba Allison Allison et al. (2006), aunque no hay consenso sobre que procedimiento es mejor para determinar el tamaño de las muestras, sí que lo hay sobre la conveniencia de realizar el análisis de la potencia, y, por supuesto, sobre el hecho de que un mayor número de replicas generalmente proporcionan una mayor potencia.
En el contexto de los microarrays, llamamos “pooling” a la combinación del RNA de diferentes casos en una única muestra.
Hay dos razones a favor de ello, una, es que, a veces, no hay suficiente ARN disponible y esta es la única forma de conseguir suficiente material para construir los arrays, otra, más controvertida, es la creencia que la variabilidad entre arrays puede reducirse por “pooling”. La justificación es que combinar muestras equivale a promediar expresiones, y como ya se conoce, el promedio es menos variable que los valores individuales. A pesar de la debilidad de este argumento, es verdad que en ciertas situaciones el pooling puede ser apropiado y, recientemente, muchos estadísticos han dedicado sus esfuerzos para tratar de responder la pregunta “to pool or not to pool” (Kerr and Churchill 2001). Por ejemplo, si la variabilidad biológica está altamente relacionada con el error en las medidas, y las muestras biológicas tienen un coste mínimo en comparación con el de los arrays, una apropiada estrategia de “pooling” puede ser claramente eficiente en costes.
De todos modos, el “pooling” no se debería usar en cualquier tipo de estudios. Si el objetivo es comparar expresiones medias (ver más adelante " comparación de clases"), puede funcionar adecuadamente, pero se debería claramente evitar cuando el objetivo del experimento es construir predictores que se basen en características individuales.
Se entiende por aleatorizar la asignación de todos los factores al azar a las unidades experimentales. Co ello se consigue disminuir el efecto de los factores no controlados por el experimentador en el diseño experimental y que podrían influir en los resultados.
Las ventajas de aleatorizar los factores no controlados son:
Hace referencia a dividir o particionar las unidades experimentales en grupos llamados bloques de modo que las observaciones realizadas en cada bloque se realicen bajo condiciones experimentales lo más parecidas posibles. A diferencia de lo que ocurre con los factores tratamiento, el experimentador no está interesado en investigar las posibles diferencias de la respuesta entre los niveles de los factores bloque.
Bloquear es una buena estrategia siempre y cuando sea posible dividir las unidades experimentales en grupos de unidades similares. La ventaja de bloquear un factor que se supone que tiene una clara influencia en la respuesta, pero en el que no se está interesado, es que convierte la variabilidad sistemática no planificada en variabilidad sistemática planificada.
En arrays de dos colores, se aplican dos condiciones experimentales a cada array. Esto permite la estimación del efecto del array, como el efecto de bloqueo. En Affymetrix u otro array de un canal, cada condición se aplica a un chip separado, imposibilitando la estimación del efecto de los microarrays, lo cual, por otra parte, se considera que tiene una relación muy pequeña en el tratamiento de los efectos debido al proceso industrial utilizado para fabricar estos chips.
Como consecuencia de lo anterior, los experimentos que usan arrays de un canal se consideran “estándar”, por lo que se les pueden aplicar los conceptos y técnicas tradicionales de diseño de experimentos .
Los arrays de dos canales presentan una situación más complicada. Por una parte, los “dos colores” no son simétricos, es decir, con la misma cantidad de material, un array hibridado con uno u otro color sea Cy5 o Cy3, emitirá señales con diferente intensidad. La forma normal de manejar este problema es dye-swapping que consiste en utilizar para una misma comparación dos arrays con dyes cambiados, es decir, si en el primer array la muestra 1 se marca con Cy3 y la muestra 2 con Cy5, con las muestras del segundo array se hace al revés (ver figura 3.3.
Figure 3.3: Representación simplificada de dos diseños
Por otra parte, el hecho de que solo se puedan aplicar dos condiciones a cada array complica el diseño, ya sea porque normalmente hay más de dos condiciones, o porque no es recomendable hibridar directamente dos muestras en un array, creando pares artificiales.
El problema de la asignación eficiente de muestras a microarrays, dado un numero de condiciones a comparar y un número fijo de arrays disponibles, ha sido estudiado de forma exhaustiva.
El diseño utilizado más comúnmente en la comunidad biológica es el diseño de referencia en el que cada condición de interés se compara con muestras tomadas de alguna referencia estándar común a todos las arrays, (ver figura 3.4, imagen (a)).
Figure 3.4: (a) diseño de referencia. (b) diseño en loop.
Los diseños de referencia permiten hacer comparaciones indirectas entre condiciones de interés. La crítica más importante a esta aproximación es que el 50 de las fuentes de hibridación se utilizan para producir la señal del grupo control, de un interés no intrínseco para los biólogos. En contraste, un diseño en loop compara dos condiciones a través de otra cadena de condiciones, por lo que elimina la necesidad de una muestra de referencia.
Los estudios realizados con microarrays, sea cual sea la tecnología en que se basan, tienen una característica común: generan grandes cantidades de datos a través de una serie de procesos, que hacen que su significado no siempre sea completamente intuitivo.
Como en todo tipo de análisis, antes de empezar a trabajar con los datos, debemos de asegurarnos de que éstos son fiables y completos y de que se encuentran en la escala apropiada para proporcionar la información que pretendemos obtener de ellos.
En el caso de los microarrays solemos distinguir dos fases previas al análisis de los datos:
Exploración y control de calidad Mediante gráficos y/o resúmenes numéricos estudiamos la estructura de los datos con el fin de decidir si parecen correctos o presentan anomalías que deben ser corregidas.
Preprocesado Incluso si son correctos, los datos “crudos” no sirven para el análisis sino que necesitan preprocesados, lo que puede interpretarse como uno o más de los procesos siguientes:
A menudo -por un cierto abuso de lenguaje- se denomina “normalización” al preprocesado descrito en las etapas anteriores.
Desde la generalización del uso de los microarrays se han desarrollado muchas formas para visualizar los datos y decidir acerca de su calidad. Algunas trabajan sobre los datos obtenidos del escáner, otras lo hacen con los datos normalizados. Algunas sirven tanto para arrays de dos colores como arrays de un color. Otras son específicas de la tecnología. Con el fin de organizar esta multitud de opciones podemos diferenciar:
La exploración y el control de calidad pueden basarse en datos de bajo o de alto nivel. En cada caso se pueden aplicar ya sean tácnicas generales de visualización de datos, como histogramas, diagramas de caja o visualización en dimensiones reducidas (PCA u otras), o bien técnicas “ad-hoc” para cada tipo de datos como el MA-Plot u otras que se discuten a continuación.
La estructura de los datos de microarrays de un color o de dos colores difiere considerablemente, tanto a nivel físico (los chips y las imágenes que de ellos se obtiene) como informático, es decir en la forma en que se representan.
Tradicionalmente los arrays de dos colores o de cDNA se realizaban de forma menos automatizada que los de un color o de Affymetrix. Esto implica que, tras obtener la imagen, el escaneado del archivo “.TIFF” resultante pudiera ser llevado a cabo mediante un software independiente como Genepix, un programa que convierte las imágenes en números y genera un archivo de información (con extensión “.gpr”) a partir del cual pueden calcularse las expresiones relativas, así como valores de calidad para cada “spot” o punto escaneado en la imagen.
Para cada imagen (o sea para cada microarray) hay un archivo “.gpr” que contiene una fila por gen y varias columnas con distintos valores, por ejemplo:
La tabla siguiente muestra lo que serían las primeras filas y columnas de un archivo “.gpr” obtenido mediante el programa Genepix.
| Gen | Señal | Fondo | Señal | Fondo | “Flag” | Otros |
|---|---|---|---|---|---|---|
| Cy5 (R) | Cy5 (Rb) | Cy3 (G) | Cy3 (Gb) | |||
| gen-1 | 3.7547 | 1.8128 | 5.0672 | 1.8496 | 1 | … |
| gen-2 | 0.8331 | 0.9175 | 1.1536 | 0.9995 | 0 | … |
| gen-3 | 9.8254 | 0.2781 | 0.6921 | 0.5430 | 1 | … |
| gen-4 | 9.1539 | 0.1918 | 3.8290 | 0.0014 | 0 | … |
| gen-5 | 4.8603 | 0.2377 | 0.5338 | 0.3335 | 0 | … |
| gen-6 | 7.8567 | 1.3941 | 1.3050 | 1.4876 | 1 | … |
| gen-7 | 2.7619 | 0.8108 | 8.4916 | 1.6518 | 0 | … |
| gen-8 | 3.6618 | 0.1918 | 0.3770 | 1.1842 | 0 | … |
| gen-9 | 4.4300 | 0.5888 | 2.8161 | 0.8707 | 1 | … |
| … | … | … | … | … | … | … |
Los valores de intensidad se convierten en una única matriz de expresión que contiene una columna por chip con los valores de intensidad relativa obtenidas por ejemplo con una “sumarización” del tipo: \[ \log\frac{R-Rb}{G-Gb}, \] y una fila por gen (mismas filas que archivos “.gpr”).
La tabla siguiente muestra lo que sería una matriz de expresión derivada de cuatro archivos “.gpr” como los de la tabla anterior.
| gen | Array-1 | Array-2 | Array-3 | Array-4 |
|---|---|---|---|---|
| gen-1 | 1.65695 | -0.10820 | 1.69515 | 8.25137 |
| gen-2 | -1.82305 | 0.11350 | 1.58807 | 0.95676 |
| gen-3 | 0.01561 | 0.47682 | -27.88036 | -1.39078 |
| gen-4 | 0.42709 | 4.29319 | -4.31366 | 30.96866 |
| gen-5 | 0.04332 | 0.24126 | -10.27675 | 0.26680 |
| gen-6 | -0.02825 | 0.68408 | 2.04163 | 0.99554 |
| gen-7 | 3.50549 | -8.04635 | 0.18286 | 0.22647 |
| gen-8 | -0.23261 | -1.08477 | 0.19582 | 1.11561 |
| gen-9 | 0.50643 | 1.52147 | 0.99092 | 0.23441 |
| … | … | … | … | … |
El resultado de escanear la imagen de un array de affymetrix es un archivo de extensión “.CEL” que, a diferencia de los arrays de dos colores, está en formato binario es decir que solo puede ser leído con programas específicos para ello.
De forma similar a los arrays de dos colores, existe un archivo “.CEL” por cada microarray.
A partir de las intensidades de los archivos “.CEL” se genera la matriz de expresión, que contiene una columna por chip con los valores de intensidad absoluta, y una fila por grupo de sondas. En el caso de arrays de affymetrix existe una gran variedad de algoritmos de “sumarización” y, según cual se utilice, se obtendrá unos u otros valores de expresión. Ahora bien, éstos serán siempre medidas absolutas, es decir independientes del resto de muestras y en una escala arbitraria.
La tabla siguiente muestra las primeras filas de una matriz de expresión sumarizada correspondiente a los primeros genes de uno de los casos resueltos.
| ID_REF | GSM188013 | GSM188014 | GSM188016 | GSM188018 |
|---|---|---|---|---|
| 1007_s_at | 15630.2 | 17048.8 | 13667.5 | 15138.8 |
| 1053_at | 3614.4 | 3563.22 | 2604.65 | 1945.71 |
| 117_at | 1032.67 | 1164.15 | 510.692 | 5061.2 |
| 121_at | 5917.8 | 6826.67 | 4562.44 | 5870.13 |
| 1255_g_at | 224.525 | 395.025 | 207.087 | 164.835 |
| … | … | … | … | … |
Los ejemplos de este capítulo se basarán en dos de los conjuntos de datos descritos en el capítulo 2.
if(!require(BiocManager)) install.packages("BiocManager")
if (!(require(CCl4))){
BiocManager::install("CCl4")
}
if (!(require(estrogen))){
BiocManager::install("estrogen")
}
if (!(require(affy))){
BiocManager::install("affy")
}
if (!(require(affyPLM))){
BiocManager::install("affyPLM")
}
Aunque, en la actualidad (año 2021) el paquete más utilizado para arrays de este tipo es el paquete oligo este ejemplo utiliza el paquete affy.
library(estrogen)
library(affy)
affyPath <- system.file("extdata", package = "estrogen")
adfAffy = read.AnnotatedDataFrame("phenoData.txt", sep="", path=affyPath)
affyTargets = pData(adfAffy)
affyTargets$filename = file.path(affyPath, row.names(affyTargets))
affyRaw <- read.affybatch(affyTargets$filename, phenoData=adfAffy)
# show(affyRaw)
actualPath <- getwd()
setwd(affyPath)
allAffyRaw <- ReadAffy()
setwd(actualPath)
El resultado de la lectura es un objeto “affyraw” de clase “affyBatch”:
class(affyRaw)
## [1] "AffyBatch"
## attr(,"package")
## [1] "affy"
print(affyRaw)
## AffyBatch object
## size of arrays=640x640 features (22 kb)
## cdf=HG_U95Av2 (12625 affyids)
## number of samples=8
## number of genes=12625
## annotation=hgu95av2
## notes=
El paquete limma es conocido como paquete para la selección de genes diferencialmente expresados pero también incluye algunas funciones para la lectura y preprocesado de arrays de dos colores.
library("limma")
library("CCl4")
dataPath = system.file("extdata", package="CCl4")
adf = read.AnnotatedDataFrame("samplesInfo.txt",
path=dataPath)
#adf
targets = pData(adf)
targets$FileName = row.names(targets)
RG <- read.maimages(targets, path=dataPath, source="genepix")
attach(RG$targets)
newNames <-paste(substr(Cy3,1,3),substr(Cy5,1,3),substr(FileName,10,12), sep="")
colnames(RG$R)<-colnames(RG$G)<-colnames(RG$Rb)<-colnames(RG$Gb)<-rownames(RG$targets)<- newNames
# show(RG)
El resultado de la lectura es un objeto de classe “RGlist”.
class(RG)
## [1] "RGList"
## attr(,"package")
## [1] "limma"
print(RG)
## An object of class "RGList"
## $R
## DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,] 47 46 40 59 56 40 45
## [2,] 1095 46 39 54 56 40 47
## [3,] 47 46 41 55 52 39 46
## [4,] 46 45 41 56 54 43 47
## [5,] 55 53 45 61 71 47 63
## DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,] 70 46 65 43 59 46 63
## [2,] 69 45 61 43 62 48 68
## [3,] 81 44 57 41 58 50 63
## [4,] 78 44 60 42 59 48 62
## [5,] 91 51 83 56 71 52 69
## DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,] 77 40 91 32
## [2,] 71 41 85 32
## [3,] 70 40 39 32
## [4,] 81 41 51 32
## [5,] 81 57 95 37
## 44285 more rows ...
##
## $G
## DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,] 77 78 45 45 87 65 128
## [2,] 4558 79 44 47 81 66 131
## [3,] 79 75 46 47 84 66 122
## [4,] 78 82 44 71 86 63 135
## [5,] 93 95 92 57 86 66 132
## DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,] 107 57 90 99 84 56 68
## [2,] 110 54 89 96 168 54 72
## [3,] 114 54 84 93 78 55 71
## [4,] 113 118 84 97 83 56 75
## [5,] 133 59 95 99 102 64 84
## DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,] 71 74 83 59
## [2,] 71 75 81 60
## [3,] 73 78 62 59
## [4,] 79 77 70 61
## [5,] 99 83 92 63
## 44285 more rows ...
##
## $Rb
## DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,] 41 39 35 47 44 36 40
## [2,] 43 39 35 46 44 36 39
## [3,] 41 40 35 48 44 35 39
## [4,] 40 40 35 49 43 35 39
## [5,] 41 40 36 48 44 35 39
## DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,] 55 39 49 37 46 42 50
## [2,] 55 37 50 37 46 41 52
## [3,] 59 39 47 37 46 41 52
## [4,] 63 39 48 36 46 41 52
## [5,] 60 39 47 37 46 41 52
## DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,] 60 36 34 30
## [2,] 59 36 35 31
## [3,] 55 35 35 30
## [4,] 58 36 34 30
## [5,] 56 36 35 30
## 44285 more rows ...
##
## $Gb
## DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,] 44 44 32 34 40 39 48
## [2,] 46 43 32 35 41 38 48
## [3,] 43 43 32 35 40 38 49
## [4,] 43 43 32 35 41 38 49
## [5,] 43 42 32 35 41 38 48
## DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,] 47 35 50 50 45 37 38
## [2,] 48 35 50 49 45 37 38
## [3,] 48 35 51 49 45 37 38
## [4,] 49 35 51 49 46 37 37
## [5,] 49 35 50 49 45 37 38
## DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,] 38 40 37 36
## [2,] 38 40 37 36
## [3,] 39 40 37 36
## [4,] 38 40 37 36
## [5,] 39 40 37 36
## 44285 more rows ...
##
## $targets
## Cy3 Cy5 RIN.Cy3 RIN.Cy5 FileName
## DMSCCl319 DMSO CCl4 9.0 9.7 251316214319_auto_479-628.gpr
## DMSCCl320 DMSO CCl4 9.0 5.0 251316214320_auto_478-629.gpr
## DMSCCl321 DMSO CCl4 9.0 2.5 251316214321_auto_410-592.gpr
## DMSCCl329 DMSO CCl4 9.0 2.5 251316214329_auto_429-673.gpr
## CClDMS330 CCl4 DMSO 9.7 9.0 251316214330_auto_457-658.gpr
## 13 more rows ...
##
## $genes
## Block Row Column ID Name
## 1 1 1 1 BrightCorner BrightCorner
## 2 1 1 2 BrightCorner BrightCorner
## 3 1 1 3 (-)3xSLv1 NegativeControl
## 4 1 1 4 A_44_P317301 AW523361
## 5 1 1 5 A_44_P386163 NM_001007719
## 44285 more rows ...
##
## $source
## [1] "genepix"
##
## $printer
## $ngrid.r
## [1] 1
##
## $ngrid.c
## [1] 1
##
## $nspot.r
## [1] 430
##
## $nspot.c
## [1] 103
Los gráficos son útiles para comprobar la calidad de los datos de microarrays, obtener información sobre cómo se deben preprocesar los datos y comprobar, finalmente, que el preprocesado se haya realizado correctamente.
Siguiendo el esquema presentado en la tabla @ref{tab:tablaDiagnosticos} se presentan a continuación los distintos gráficos utilizados con una breve descripción de lo que representa cada uno y como interpretarlos adecuadamente. El código para generarlos se presenta al final del capítulo como un apéndice.
Estos gráficos permite hacerse una idea de si las distribuciones de los distintos arrays son similares en forma y posición. La figura ?? muestra los histogramas correspondientes a los 9 arrays del conjunto de datos .
affySampleNames <- rownames(pData(allAffyRaw))
affyColores <- c(1,2,2,3,3,4,4,8,8)
affyLineas <- c(1,2,2,2,2,3,3,3,3)
hist(allAffyRaw, main="Signal distribution", col=affyColores, lty=affyLineas)
legend (x="topright", legend=affySampleNames , col=affyColores, lty=affyLineas, cex=0.7)
Como los histogramas los diagramas de caja –basados en los distintos cuantiles de las valores– dan una idea de la distribución de las intensidades. La figura ?? muestra los diagramas de caja correspondientes a los 9 arrays del conjunto de datos .
boxplot(allAffyRaw, main="Signal distribution", col=affyColores, las=2)
El análisis de componentes principales puede servir para detectar si las muestras se agrupan de forma “natural” es decir con otras muestras provenientes del mismo grupo o si no hay correspondencia clara entre grupos experimentales y proximidad en este gráfico. Cuando esto sucede no significa necesariamente que haya un problema pero puede ser indicativo de efectos técnicos -como el conocido efecto “batch”- que podría ser necesario corregir.
plotPCA <- function ( X, labels=NULL, colors=NULL, dataDesc="", scale=FALSE)
{
pcX<-prcomp(t(X), scale=scale) # o prcomp(t(X))
loads<- round(pcX$sdev^2/sum(pcX$sdev^2)*100,1)
xlab<-c(paste("PC1",loads[1],"%"))
ylab<-c(paste("PC2",loads[2],"%"))
if (is.null(colors)) colors=1
plot(pcX$x[,1:2],xlab=xlab,ylab=ylab, col=colors,
xlim=c(min(pcX$x[,1])-10, max(pcX$x[,1])+10),
ylim=c(min(pcX$x[,2])-10, max(pcX$x[,2])+10),
)
text(pcX$x[,1],pcX$x[,2], labels, pos=3, cex=0.8)
title(paste("Plot of first 2 PCs for expressions in", dataDesc, sep=" "), cex=0.8)
}
if (!file.exists("allAffyNorm")){
allAffyNorm<- rma(allAffyRaw)
affyNorm <- rma(affyRaw)
save(allAffyNorm, affyNorm, file="affyNorm.Rda")
}else{
load(file="affyNorm.Rda")
}
La figura ?? muestra dos diagramas de compnentes principales realizados a partir de los datos normalizados del conjunto de datos . El gráfico de la parte superior que incluye el array defectuoso ilustra que la principal fuente de variabilidad es la diferencia de este array con el resto. Cuando se repite el análisis omitiendo esta muestra puede verse como la principal fuente de variación (eje X) se asocia con el tiempo de exposición (alto a la derecha, bajo (10h) a la izquierda, mientras que la segunda fuente de variación se asocia con la exposición a los estrógenos (alto arriba, bajo abajo).
El gráfico de la parte inferior en el que este array se ha eliminado antes de calcular las componentes principales muestra dos agrupaciones claras de arrays “derecha e izquierda” y “arriba y abajo” de la figura. El hecho de que los grupos no se correspondan con los tratamientos sugiere que puede haber alguna fuente adicional de variación que ha de tenerse en cuenta.
opt <- par(mfrow=c(2,1))
plotPCA(exprs(allAffyNorm), labels=affySampleNames, dataDesc="PCA for all arrays\nincludes defective sample")
plotPCA(exprs(affyNorm), labels=colnames(exprs(affyNorm)), dataDesc="PCA for all arrays")
par(opt)
Otra forma de ver si las muestras se agrupan según los grupos experimentales, o mediante otros criterios es usando un cluster jerárquico que realiza una agrupación básica de las muestras por grado de similaridad según la distancia que se utilice.
Como en el caso de las componentes principales si las muestras se agrupan según las condiciones experimentales es una buena señal pero si no es así puede deberse a la presencia de otra fuente de variación o bien al hecho de que se trata de un gráfico basado en todo los datos y las condiciones experimentales pueden haber afectado un pequeño número de genes.
La figura 4.1 muestra como se agrupan los datos del conjunto en base a un cluster jerárquico. Como en el caso de las componentes principales tras eliminar el array defectuoso las muestras se separan, primero por el tiempo de exposición y luego por niveles de estrógeno suministrado.
clust.euclid.average <- hclust(dist(t(exprs(affyNorm))),method="average")
plot(clust.euclid.average, labels=colnames(exprs(affyNorm)), main="Hierarchical clustering of samples", hang=-1, cex=0.7)
Figure 4.1: Un cluster jerárquico sirve para determinar si las muestras se agrupan de forma natural según los grupos experimentales o si lo hacen de otra forma
El diagnóstico de arrays de dos canales se basa principalmente en la imagen y en diferentes tipos de gráficos.
La normalización discutida en este mismo capítulo es un punto clave en el proceso de análisis de microarrays y se ha dedicado un gran esfuerzo a desarrollar y probar diferentes métodos (Quackenbush (2002),Yang et al. (2002)). Una razón para ello es que hay diferentes artefactos técnicos que deben ser corregidos para poder ser utilizados, y no cualquier método puede funcionar con todos ellos.
En general, los métodos de normalización se basan en el siguiente principio: .
Esto da una idea de como debería ser un gráfico de intensidades. Por ejemplo, si no hubiese artefactos técnicos, en arrays de dos canales, una gráfica de dispersión de intensidad del rojo frente al verde debería dejar la mayor parte de los puntos alrededor de una diagonal. Cualquier desviación de esta situación debería ser atribuible a razones técnicas, no biológicas, y por tanto, debería ser eliminada. Esto ha conducido a un método de normalización muy popular consistente en estimar la transformación a aplicar, como una función de las intensidades utilizando el método en la representación transformada de la gráfica de dispersión conocida como el .
La figura 4.2 muestra, en su parte izquierda (a) un gráfico de dispersión del canal rojo frente al verde en un array de dos colores. El hecho de que los datos no estén centrados alrededor de la diagonal sugiere la necesidad de normalización.
Una representación muy popular que ayuda a visualizar mejor esta asimetría es lo que se conoce como “”, que aparece en la parte derecha (b) de la figura4.2. Geométricamente representa una rotación del gráfico de dispersión, en la que el significado de los nuevos ejes es:
Figure 4.2: (a) Gráfico de un canal frente al otro (b) MA-plot (intensidad frente log-ratio)
La imagen del chip (véase la figura ??, izquierda) ofrece una visión rápida de la calidad del array, proporcionando información acerca del balance del color, la uniformidad en la hibridación y en los , de si el background es mayor del normal y dela existencia de artefactos como el polvo o pequeñas marcas (rasguños).
Estos gráficos (véase la figura 4.3, derecha) son útiles para detectar posibles anormalidades o un background excesivamente alto .
Figure 4.3: La relación señal ruido sirve para detectar posibles anormalidades o un background excesivamente alto como medida de calidad
Un gráfico muy utilizado es el diagrama de cajas o “boxplot” múltiple con una caja por cada chip. Del alineamiento (o falta de él) y la semejanza (o disparidad) entre las cajas, se deduce si hace falta, o no, normalizar entre arrays.
En el caso de arrays de dos colores pueden utilizarse diagrams de cajas “dentro de arrays” (entre distintos sectores del mismo chip) y “entre arrays”.
Los arrays de affymetrix contienen millones de sondas por lo que no pueden examinarse a simple vista. A pesar de ello hay diversas formas de obtener una imagen que, en caso de presentar irregularidades pueden indicar algún tipo de problemas como burbujas, arañazos, etc. La figura @ref(c06plotAffy} muestra algunas de las imágenes
La imagen del array de Affymetrix sólo es útil para evidenciar grandes problemas como burbujas, arañazos, etc. En este ejemplo los dos arrays de la izquierda se considerarían aceptables y los de la derecha defectuosos.
Figure 4.4: Imágenes de cuatro microarrays de Affymetrix
En los chips de dos colores el MA–plot se utliza para comparar los dos canales en cada array (rojo y verde). En cambio, en los chips de Affymetrics, en que sólo hay un canal en cada array, la única forma de definir M (el log ratio) es a partir de la comparación entre pares de de valores, ya sea los arrays dos a dos o bien cada array respecto un valor de referencia que puede ser la mediana, punto a punto, de todos los arrays (véase por ejemplo la figure ).
Figure 4.5: En los chips de Affymetrix la única forma de definir M (el log ratio) es comparar entre diferentes arrays
\(M=\log_2(I_1) - \log_2(I_2)\): log ratio \(A=\displaystyle \frac{1}{2}(\log_2 (I_1)+\log_2(I_2))\): log de intensidades Donde \(I_1\) es la intensidad del array de estudio, e \(I_2\) es la intensidad media de arrays. Por lo general, se espera que la distribución en el gráfico se concentre a lo largo del eje M = 0.
Los modelos de bajo nivel (“Probe-Level-Models” o PLM) ajustan a los valores de intensidad –a nivel de sondas, no de valores totalizados de gen– un modelo explicativo. Los valores estimados por este modelo se comparan con los valores reales y se obtienen los errores o “residuos” del ajuste. El análisis de dichos residuos procede de forma similar a lo que se realiza al analizar un modelo de regresión: Si los errores no presentan ningún patrón especial supondremos que el modelo se ajusta relativamente bien. Si, en cambio, observamos desviaciones de esta presunta aleatoriedad querrá decir que el modelo no explica bien las observaciones, lo cual se atribuirá a la existencia de algún problema con los datos.
Con los valores ajustados del modelo se calculan dos medidas:
stopifnot(require(affyPLM))
Pset<- fitPLM(allAffyRaw)
La figura ?? muestra los gráficos RLE y NUSE para el conjunto de datos estrogen. En ambos gráficos puede verse como el array defectuoso “” queda claramente diferenciado del resto.
opt<- par(mfrow=c(2,1))
RLE(Pset, main = "Relative Log Expression (RLE)",
names=rownames(pData(allAffyRaw)), las=2, cex.axis=0.6)
NUSE(Pset, main = "Normalized Unscaled Standard Errors (NUSE)", las=2,
names=rownames(pData(allAffyRaw)), las=2, cex.axis=0.6)
(#fig:plotPLM )Graficos de diagnóstico calculados a nivel de sondas PLM
par(opt)
La palabra describe las técnicas utilizadas para transformar adecuadamente los datos antes de que sean analizados. El objetivo es corregir diferencias sistemáticas entre muestras, en la misma o entre imágenes, lo que no representa una verdadera variación entre las muestras biológicas.
Estas diferencias sistemáticas pueden deberse, entre otras, a:
A veces puede ser difícil detectar estos problemas , aunque existen algunas formas de saber si es necesaria realizar una normalización. Aqui destacamos dos posibilidades:
Este método esta basado en un ajuste global, es decir en modificar todos los valores una cantidad , estimada de acuerdo a algún criterio. \[\begin{equation} \log_2 R/G \rightarrow \log_2 R/G-c=\log_2 R/(Gk) \end{equation}\] opciones para \(k\) o \(c= \log_2k\) son
\(c\)= mediana o media de log ratio para un conjunto concreto de genes o genes control o genes housekeeping.
La intensidad total de la normalización
\(k=\sum R_i/\sum G_i\)
En este caso se realiza una modificación específica para cada valor. Esta modificación se obtiene como una función de la intensidad total del gen (\(c=c(A)\)). \[\begin{equation} \log_2 R/G \rightarrow \log_2 R/G-c(A)=\log_2 R/(Gk(A)) \end{equation}\]
Una posible estimación de esta función puede hacerse utilizando la función (LOcally WEighted Scatterplot Smoothing).
En los arrays de Affymetrix, como en todos los tipos de microarrays, tras escanear la imagen se obtiene una serie de valores de intensidad de cada elemento del chip. En el caso de estos arrays sabemos que cada valor no corresponde a la expresión de un gen:
El proceso que convierte las señales individuales en valores de expresión normalizados para cada gen consta de tres etapas:
A menudo los tres pasos se denominan genérica -y erróneamente- “normalización”.
A diferencia de los chips de ADNc, aquí las medidas de expresión son absolutas (no se compara una condición contra otra) dado que cada chip se hibrida con un única muestra.
Hay muchos métodos para estimar la expresión (más de 30 publicados). Cada método contempla de forma explícita o implícita las tres formas de preprocesado: corrección del fondo, normalización y resumen.
Los principales métodos que consideraremos son:
Es el primer método introducida por Affymetrix. La corrección del fondo se realiza restando el “perfect match” del “mismatch” \[\begin{equation} E_j=PM_j-MM_j \end{equation}\] La normalización se realiza de forma global haciendo transformaciones de forma que la media de todo el chip sea la misma y la sumarización se basa en calcular el promedio de las diferencias absolutas ignorando los pares que se desvían más de \(3\sigma\) de \(\mu\). \[\begin{equation} Dif. Media=\frac{1}{|A|}\sum_{j \in A}(PM_j-MM_j) \end{equation}\] Los problemas que presenta estemétodo son:
Los problemas que presentaba el método M.A.S. 4.0 llevaron a sustituirlo por otra variante, el M.A.S 5.0, llamado así por venir implementado en el software de affymetrix llamado “MicroArray Suite 5.0”.
Este método utiliza un estadístico robusto, , para corregir y ponderar el fondo y calcular (estimar) la señal. El biweight de Tukey \(T_{bi}\) pondera los valores por su distancia a la mediana, es decir, mide la tendencia central pero realiza un ajuste de outliers.
La lógica de este método reside en pensar que el valor de MM no siempre tiene sentido, (p.ej si MM \(>\) PM). Dado que esto sucede en ocasiones se realiza el cambio siguiente:
Este método no tan solo corrige el background sino que también permite normalizar y sumarizar. Para ello se introduce el “Mismatch Idealizado” (IM) que permite corregir la intensidad de las pruebas individuales. Este método también ha sido muy criticado:
Para compensar algunas deficiencias de los primeros métodos de resumen y normalización de arrays de Affymetrix, Irizarry y sus colegas introdujeron en 2003 (~) un método basado en la modelización de las intensidades de las sondas que, en vez de basarse en las distintas sondas de un gen dentro de un mismo array se basa en los distintos valores de la misma sonda entre todos los arrays disponibles,
Esquemáticamente los pasos que realiza este método son:
Como resultado final de todos los pasos anteriores se obtiene la matriz con los datos sumarizados y normalizados. A pesar de no estar exento de críticas como la que afirma que este procedimiento “compacta” los valores reduciendo su variabilidad natural, este método se ha convertido en el estándar “de facto” actualmente por muchos usuarios de Bioconductor.
Figure 4.6: El método RMA incluye una normalización por cuantiles como la representada esquemáticamente en esta figura
El filtraje no específico es recomendable para eliminar el ruido de fondo y limitar los ajustes posteriores a los necesarios. Los principales procesos de filtrado son:
El objetivo del filtraje es eliminar aquellos spots cuyas imágenes o señales sean erróneas por diferentes motivos, disminuyendo el ruido de fondo. Aunque existe controversia a su uso, prefiriendo el no filtrado a eliminar de forma no intencionado spots informativos.
El motivo más habitual para el que se suelen utilizar microarrays es la búsqueda de genes cuya expresión cambia entre dos o más condiciones experimentales, por ejemplo a consecuencia de un tratamiento, una enfermedad u otras causas (distintos tiempos, distintas lineas celulares, …).
El problema consiste en identificar estos genes y suele denominarse o bien comparación de clases.
El problema de seleccionar genes diferencialmente expresados se traduce de manera casi inmediata al problema estadístico de comparar variables y, en años recientes, se han desarrollado un gran número de métodos estadísticos para resolverlo. La mayoría son extensiones de los métodos estadísticos clásicos, pruebas-t o análisis de la varianza, adaptados en uno u otro sentido para tener en cuenta las peculiaridades de los microarrays.
Aunque el problema de la selección de genes diferencialmente expresados puede relacionarse directamente con la realización de pruebas estadísticas, en el caso de los microarrays, el hecho de que haya dos tecnologías que miden la expresión de dos formas distintas hace que se deba diferenciar la metodología a emplear en cada caso. Los arrays de dos colores combinan dos muestras en un chip y generan una medida de expresión relativa. Esto hace que para comparar dos muestras de un mismo individuo sean la opción naturalmente más apropiada
En el caso de querer comparar muestras independientes de diferentes individuos los arrays de un color son la mejor opción. Evidentemente, lo más común será disponer de una sola técnica y tener que adaptar los análiaisis estadísticos a la misma.
Vamos a plantear un posible esquema de trabajo, para situar la mejor opción en cada caso:
Recordemos que una vez se han hecho los experimentos con microarrays y obtenida la señal, los datos que se disponen son los logaritmos del valor detectado por el escaner.Esto hace que algunas operaciones que se realicen tengan en cuenta esta característica. Segun si la comparación a realizar se llevará a cabo con datos independientes (2 muestras, casos 1 y 2) o con datos dependientes (muestras apareadas, casos 3, 4) algunas medidas naturales o razonables para la comparación de expresiones son las siguientes:
Consideremos la tabla siguiente que representa una matriz de expresión simplificada que contiene las expresiones relativas (por ejemplo entre tejido tumoral y sano del mismo individuo) de 5 genes en 6 muestras.
Podemos calcular las medidas descritas para el caso de una muestra para decidir si un gen está expresado o no lo está. Se discutirá más adelante como precisar esto pero de momento nos quedaremos con la idea de que si la medida escogida es (cercana a) cero el gen no está diferencialmente expresado y si es mayor o menor que cero si que lo está. Nos referimos a “cero” porque estamos hablando de logaritmos de razones: si la expresión es la misma en ambas condiciones el cociente es uno y su logaritmo es cero.
Vale la pena insistir en el concepto de expresión diferencial: no nos preocupa cual es la expresión del gen en una u otra muestra sinó si son distintas.
La tabla anterior sugiere que podría considerarse el gen A está diferencialmente expresado (promedio y t–test altos) mientras que el gen B o el D no lo están (promedio y test-t próximos a cero). Los genes C y D pueden llevar a conclusiones contradictoria según nos basemos en el promedio o el test t.
Si se observan los valores del test t del gen C se concluye que el gen no aparenta estar diferencialmente expresado. Si en cambio se observa su promedio parece que si que lo esté.
En el gen E pasa exctamente lo opuesto. La explicación de estas aparentes contradicciones se halla en el error estándar. En el gen C es muy elevado, debido a que el valor (20) es probablemente un “outlier”. En el gen D el error estándar es muy bajo por lo que, al encontrarse en el denominador del t-test aumenta artificialmente su valor.
Vamos a plantear la forma como se aborda este problema en el estudio de datos de microarrays. Dadas las características propias de este tipo de datos se consideran otras formas de estimar el error estándar que no sean tan sensibles a valores extremos o muy bajos.
Consideremos el gen \(g\). Si llamamos:
Podemos considerar dos variantes para el test-t:
| Test |
| Test-t Global |
| Test-t específico |
\end{center}
En la práctica muchos métodos de selección de genes diferencialmente expresados han acabado buscando un compromiso entre ambas aproximaciones para lo que proponen o derivan fórmulas que de alguna forma ponderan o combinan dos estimaciones del error estándar: una basada en todos los genes y otra específica de cada gen. La tabla siguiente ilustra como algunos de los métodos más utilizados en la bibliografía incorporan esta idea.
Al hecho de incluir un coeficiente que tenga en cuenta la variabilidad de todos los genes en el array para estimar el error estándar de cada gen se le denomina moderación de la varianza (“variance shrinkage”) y es una de las aproximaciones en que existe cierto consenso (Allison et al. (2006)) acerca de que sirven para mejorar la selección de genes diferencialmente expresados.
Como ejemplo utilizaremos el conjunto de datos celltypes y supondremos que disponemos ya de los datos normalizados y filtrados almacenados en un objeto expressionSet.
Este proceso es sencillo de hacer siguiendo los pasos del capiítulo anterior. Para simplificar el ejemplo una vez cargados y normalizados filtraremos los datos usando funcion nSFilter de la librería genefilter con sus opciones por defecto, lo que nos dejara un total de 10254 sondas en vez de las 45101 originales.
library(affy)
library(genefilter)
affyPath<- "datos/celltypes/celfiles"
cellTypesTargets<- read.AnnotatedDataFrame("targets.txt", sep="\t", path=affyPath)
cellTypesTargets$filename = file.path(affyPath, row.names(cellTypesTargets))
cellTypesRaw <- read.affybatch(cellTypesTargets$filename, phenoData=cellTypesTargets)
eset_rma <- rma(cellTypesRaw)
Filtered <- nsFilter(eset_rma)
eset_rma_filtered <- Filtered[["eset"]]
save(eset_rma, eset_rma_filtered, file="datos/celltypes/celltypes-normalized.rma.Rda")
Si consideramos que el campo treat del objeto contiene la información de los grupos a comparar (ratones estimulados con LPS frente a los que no han sido tratados), podemos realizar un primer análisis utilizando el test-t. Para ello utilizaremos el la función rowttests del paquete genefilter que realiza un test \(t\) sobre cada una de las filas (genes) de una matriz de expresión.
stopifnot(require(Biobase))
load (file="./datos/celltypes/celltypes-normalized.rma.Rda")
my.eset <- eset_rma_filtered
grupo_1 <- as.factor(pData(my.eset)$treat)
stopifnot(require(genefilter))
teststat <-rowttests(my.eset, "treat")
print(teststat[1:5,])
## statistic dm p.value
## 1424378_at -10.6739944 -1.1113660 8.715342e-07
## 1417675_a_at 20.7390631 0.9451576 1.504712e-09
## 1436530_at 9.2313983 1.4943922 3.291690e-06
## 1428603_at -3.4882096 -0.4013432 5.840466e-03
## 1458888_at 0.7083721 0.2343463 4.948949e-01
Cuanto mayor sea el valor absoluto del estadístico \(t\) mayor es la probabilidad de que el gen esté diferencialmente expresado.
Como hemos visto en la sección anterior dos medidas naturales para la selección de genes son el promedio de “log-ratios”, o la diferencia de promedios en el caso de muestras independientes, o el valor del estadístico de test (\(t\)-test) de una o dos muestras según si se trata de muestras apareadas o independientes respectivamente. Los primeros estudios de microarrays eran muy costosos y se hacían con pocas o incluso ninguna réplica por condición experimental. En estas situaciones la única forma fiable de detectar una diferencia de expresión era a traves del “log-ratio” o su diferencias.
Rápidamente se puso en evidencia que para poder obtener los genes que estaban realmente diferencialmente expresados era preciso disponer de un soporte estadístico que permitiera tener en cuenta la variación aleatoria existente entre muestras.
En la práctica esto se reduce a afirmar que si, además de la diferencia de expresión entre las condiciones experimentales, se lleva a cabo un test estadístico dispondremos de una medida objetiva, el p–valor que nos servirá para decidir qué genes se declaran diferencialmente expresados, a saber, aquellos en los que el p–valor del test sea inferior a un cierto umbral como 0.05 o 0.01.
Como es sabido, un test estadístico procede decidiendo rechazar la hipótesis nula si el p–valor es más pequeño que el nivel de significación del test.
Siguiendo con el ejemplo anterior podemos ordenar los resultados de los tests en base a los p–valores:
ranked <-teststat[order(teststat$p.value),]
print(ranked[1:5,])
## statistic dm p.value
## 1449383_at -48.61014 -2.044928 3.276616e-13
## 1451421_a_at -40.72173 -1.445182 1.909283e-12
## 1415929_at -39.77036 -0.812879 2.415238e-12
## 1450826_a_at 37.62239 5.147992 4.193941e-12
## 1448303_at -31.92365 -2.283765 2.140612e-11
Ahora podríamos seleccionar, por ejemplo los genes cuyo p–valor fuera inferior a 0.01
selectedTeststat <- ranked[ranked$p.value < 0.01,]
Esto deja un total de {r} dim(selectedTeststat)[1] con un p–valor inferior a 0.01
Si se opta por computar los valores de significación (p–valores) de los genes, resulta interesante comparar el tamaño del cambio del nivel de significación estadístico. El “volcano plot” es una representación gráfica que permite ordenar los genes a lo largo de dos dimensiones, la biológica, representada por el “fold change” y la estadística representada por el logaritmo negativo del p–valor.
En la escala horizontal se representa el cambio entre los dos grupos (en escala logarítmica, de manera que la regulación positiva o negativa se representa de forma simétrica). En la escala vertical se representa el p–valor del test en una escala logarítmica negativa, de forma que los p–valores más pequeños aparecen mayores.
Así pues puede considerarse que el primer eje indica impacto biológico del cambio (a más efecto biológico mayor “fold-change”) y el segundo muestra la evidencia estadística, o la fiabilidad del cambio.
Como veremos más adelante en esta sección, el paquete limma permite realizar volcano-plots de forma sencilla a partir de los resultados de un anàlisis, pero para hacer un volcano plot tan sólo se necesita una lista de p-valores y los efectos (“fold-change”) asociados.
La figura (???)(fig:volcano0) muestra un “volcano-plot” para este ejemplo de los “celltypes” para la comparación “LPS” frente a “Medium”.
FC <- teststat$statistic
pVal<-teststat$p.value
Y <- -log(pVal)
plot(Y~FC)
También podemos adaptar alguna función “ad-hoc” como la siguiente, tomada de :
library(ggrepel)
# plot adding up all layers we have seen so far
volcanoP<- function (de,log2FoldChange, pvalue,
diffexpressed=NULL, col=NULL, delabel=NULL){
ggplot(data=de, aes(x=log2FoldChange, y=-log10(pvalue), col=diffexpressed, label=delabel)) +
geom_point() +
theme_minimal() +
geom_text_repel() +
scale_color_manual(values=c("blue", "black", "red")) +
geom_vline(xintercept=c(-0.6, 0.6), col="red") +
geom_hline(yintercept=-log10(0.05), col="red")
}
diffexpressed <- ifelse(abs(teststat$statistic)>10, TRUE, FALSE)
label <- rep(NA, length(diffexpressed))
label[diffexpressed] <- rownames(teststat)[diffexpressed]
volcanoP (de=teststat, log2FoldChange=teststat$statistic, pvalue=teststat$p.value,
diffexpressed = diffexpressed, delabel = label)
Tal como se ha indicado más arriba, para realizar un suele controlarse la probabilidad de error de tipo I (de falsos positivos) y buscar, de entre los tests candidatos, aquellos que tengan una menor probabilidad de error de tipo II, o equivalentemente una mayor potencia. A partir de este planteamiento existe, en el contexto estadístico estándar, multitud de formas de determinar cual debe ser el tamaño muestral necesario para obtener una potencia dada fijados el tamaño de efecto (“fold-change”) y la probabilidad de error de tipo I.
En el caso de los microarrays se han desarrollado diversas fórmulas para realizar cálculos de este tipo, pero la compleja estructura de los datos microarrays hace que sean relativamente discutibles.
Simon (Simon et al. (2003)) sugiere la fórmula siguiente que es una generalización de las fórmulas clásicas para problemas de dos muestras:
El tamaño total requerido para detectar genes diferencialmente expresados en al menos una diferencia \(\delta\) con una probabilidad de error de tipo I (FP), \(\alpha\) y una probabilidad de error de tipo II (FN) \(1-\beta\) se calcula: \[ n=\frac{4(z_{\alpha/2}+z_\beta)^2}{(\delta/\sigma)^2}, \] donde \(z_\alpha\) y \(z_\beta\) son los percentiles 100\(\alpha/2\) y 100\(\beta\) de la distribución Normal \(N(0,1)\) y \(\sigma\) es la desviación estandar de un gen dentro de una clase (de un grupo). Obviamente \(\sigma\) es siempre desconocida por lo que, sin una muestra piloto con que estimarla el calculo e más imaginativo que realista.
Además de esto, el número de arrays usualmente recomendado queda lejos de la cantidad asequible para la mayor parte de los experimentos ((???),Tibshirani (2006)). Lo que muchos usuarios hacen a la práctica, es buscar un equilibrio entre los costes y la reproducibilidad y, de hecho, tienden a usar una cantidad fija de arrays tal como 3 o 5 sin consideraciones adicionales.
Por ejemplo si ponemos \(\alpha=0.001\), \(1-\beta=0.95\), \(\delta=1\) y estimamos \(\sigma\) entre todas las muestras el número de réplicas biológicas que necesitaremos será de 35.8.
El análisis de microarrays se realiza en base gen a gen e involucra múltiples tests, miles probablemente. Esto significa que, a medida que crece el número de genes, la probabilidad de declarar erroneamente al menos un gen diferencialmente expresado va en aumento, y si no se realiza algún tipo de ajuste el número de falsos positivos será tanto más alto cuantos más genes estemos analizando.
Hay muchas formas de intentar controlar estas probabilidades de error y puede verse un excelente revisión en Dudoit Dudoit, Shaffer, and Boldrick (2003)).
De forma simplificada consideramos las dos aproximaciones más utilizadas.
Una posibilidad es mirar de controlar la probabilidad de obtener al menos un falso positivo o “Family-wise-error-rate (FWER)”. El más popular de estos métodos de control es la corrección de Bonferroni, consistente en multiplicar el p–valor por el número de tests realizados. La misma Dudoit y muchos otros autores han desarrollado variantes de los métodos clásicos de ajuste FWER usando por ejemplo tests de permutaciones.
El criterio FWER es quizás demasiado restrictivo dado que el control de los falsos positivos implica un considerable incremento de falsos negativos. En la práctica, sin embargo, muchos biólogos parecen estar dispuestos a aceptar que se produzcan algunos errores, siempre y cuando esto permita realizar descubrimientos. Por ejemplo un investigador debe considerar aceptable cierta pequeña proporción de errores (digamos del 10 al 20%) entre sus descubrimientos. En este caso, el investigador está expresando interés en controlar el porcentaje de falsos descubrimientos (FDR), es decir lo que es la proporción de falsos positivos sobre el total de genes inicialmente identificados como expresados diferencialmente. A diferencia del nivel de significación que queda determinado antes de examinar los datos, la FDR es una medida de confianza a posteriori ya que emplea información disponible en los datos para estimar las proporciones de falsos positivos que han tenido lugar. Si se obtiene una lista de los genes expresados diferencialmente en los que el FDR se controla hasta, digamos, el 20%, cabe esperar que el 20% de estos genes representen falsos positivos. Lo cual supone un enfoque menos restrictivo que controlar el FWER.
La decisión de controlar el FDR o el FWER depende de los objetivos del experimento. Si, por ejemplo, el objetivo es la es razonable permitir cierta cantidad de falsos positivos y es preferible seleccionar FDR. Si por el contrario se trabaja con una lista de un tamaño menor al deseado para verificar la expresión de ciertos genes específicos, entonces el FWER es el criterio apropiado.
El ejemplo siguiente muestra como se realiza el ajuste de p–valores usando el paquete multtest de Bioconductor para ajustar por los métodos de Bonferroni (FWER), Benjamini & Hochberg (FDR) o Benjamini & Yekutieli (BY, FDR).
stopifnot(require(multtest))
procs <- c("Bonferroni","BH", "BY")
adjPvalues <- mt.rawp2adjp(teststat$p.value, procs)
names(adjPvalues)
## [1] "adjp" "index" "h0.ABH" "h0.TSBH"
ranked.adjusted<-cbind(ranked, adjPvalues$adjp)
head(ranked.adjusted)
## statistic dm p.value rawp Bonferroni
## 1449383_at -48.61014 -2.044928 3.276616e-13 3.276616e-13 3.359842e-09
## 1451421_a_at -40.72173 -1.445182 1.909283e-12 1.909283e-12 1.957778e-08
## 1415929_at -39.77036 -0.812879 2.415238e-12 2.415238e-12 2.476585e-08
## 1450826_a_at 37.62239 5.147992 4.193941e-12 4.193941e-12 4.300467e-08
## 1448303_at -31.92365 -2.283765 2.140612e-11 2.140612e-11 2.194983e-07
## 1418645_at -28.74611 -4.126867 6.044555e-11 6.044555e-11 6.198087e-07
## BH BY
## 1449383_at 3.359842e-09 3.296908e-08
## 1451421_a_at 8.255283e-09 8.100651e-08
## 1415929_at 8.255283e-09 8.100651e-08
## 1450826_a_at 1.075117e-08 1.054978e-07
## 1448303_at 4.389966e-08 4.307737e-07
## 1418645_at 1.033014e-07 1.013665e-06
Si seleccionamos los genes en base a su p–valor ajustado por ejemplo por le método de Banjamini y Yekutieli se obtienen los siguientes genes
selectedAdjusted<-ranked.adjusted[ranked.adjusted$BY<0.001,]
nrow(selectedAdjusted)
## [1] 420
En la sección anterior se ha discutido el uso del test \(t\) y sus extensiones para la selección de genes diferencialmente expresados en situaciones relativamente sencillas, es decir cuando sólo hay dos grupos.
En muchos estudios el número de grupos a considerar es más de dos y las fuentes de variabilidad pueden ser más de una, por ejemplo una puede ser el tratamiento pero otra puede ser la edad de los individuos o cualquier otra condición fijada por el investigador o derivadad de la heterogeneidad de las muestras.
En estas situaciones una aproximación razonable en problemas con una variable respuesta es el análisis de la varianza, discutido en el capítulo here. En esta sección se presenta una aproximación equivalente que de forma general se denomina el modelo lineal. Este método -que engloba el análisis d la varianza y la regresión- es uno de los más usados en estadística y ha sido popularizado en el campo de microarrays gracias a los trabajos de Gordon Smyth quien ha creado el paquete limma que se ha convertido en la herramiente más utilizada para el análisis de microarrays.
El modelo lineal (ver por ejemplo Faraway, (???)) es un marco general para la modelización y el análisis de datos estadística.
EL método consiste en asumir una relación lineal entre los valores observados de una variable respuesta y las condiciones experimentales. A partir de aquí se obtienen estimadores para los parámetros del modelo y de sus errores estándar, y (con algunas condiciones extra) es posible hacer inferencia acerca del experimento.
La aplicación de modelos lineales puede ser visto como un proceso secuencial, con los siguientes pasos:
El esquema anterior se puede aplicar a casi cualquier tipo de situación experimental. En la sección siguiente se presentan un par de ellas.
Swirl es una mutación puntual que provoca defectos en la organización del embrión en desarrollo a lo largo de su eje dorsal-ventral. Como resultado, algunos tipos celulares se reducen y otros se expanden. Un objetivo de este trabajo fue identificar los genes con expresión alterada en el mutante Swirl en comparación con “wild zebrafish”.
| 1 |
| 2 |
| 3 |
| 4 |
\end{center}
Los plásmidos IncHI codifican genes de resistencia múltiple a los antibióticos en S. enterica.
El plásmido R27 de la cepa salvaje es termosensible al transferirse.
Algunos fenotipos mutantes relacionados con hha y hns cromosómicos participan en diferentes procesos metabólicos de interés en la conjugación termoregulada.
El objetivo del experimento es encontrar qué genes se expresean diferencialmente en tres tipos de mutantes diferentes, \(M_1\), \(M_2\) y \(M_3\).
Este experimento debe ser planteado de forma diferente según el tipo de arrays (uno o dos colores) y qué comparaciones son las de mayor interés.
" width=“50%”>
" width=“50%”>
" width=“50%”>
Modelo, \(\mathbf{y=X\alpha+\varepsilon}\), y contrastes \(\mathbf{C'\beta}\)
\[ \left( \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \end{array} \right) = \underbrace{\left( \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ -1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & -1 \end{array} \right)}_{\text{Matriz de Diseño}, \mathbf{X}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \end{array} \right) + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6 \end{array} \right) \] \[ \left( \begin{array}{r} \beta_1 \\ \beta_2 \\ \beta_3 \end{array} \right) = \underbrace{\left( \begin{array}{rrr} 1 & -1 & 0 \\ 1 & 0 & -1 \\ 0 & 1 & -1 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \end{array} \right). \]
Modelo, \(\mathbf{y=X\alpha+\varepsilon}\), y contrastes \(\mathbf{C'\beta}\)
Parámetros del modelo: \[ \alpha_1= \mathbf{E}\left (log\frac{M_1}{M_2}\right ),\ \alpha_2= \mathbf{E}\left (log\frac{M_2}{M_3}\right ). \] \(\alpha_3\) no es necesaria: \(\log \left( \frac{M_1}{M_3}\right )=\log\left(\frac{M_1}{M_2}\right )-\log\left(\frac{M_2}{M_3}\right)\).
Contrastes: Algunas de las comparaciones deseadas son precisamente los parámetros. \[\begin{array}{ccc} \beta_1 &=& \alpha_1,\\ \beta_2 &=& \alpha_2,\\ \beta_3 & =& \alpha_1+\alpha_2. \end{array}\]
\[ \left( \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \end{array} \right) = \underbrace{\left( \begin{array}{rr} 1 & 0 \\ 0 & 1 \\ 1 & -1 \\ -1 & 0 \\ 0 & -1 \\ -1 & 1 \end{array} \right)}_{\text{Matriz de Diseño}, \mathbf{X}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \end{array} \right) + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6 \end{array} \right) \]
\[ \left( \begin{array}{r} \beta_1 \\ \beta_2 \\ \beta_3 \end{array} \right) = \underbrace{\left( \begin{array}{rrr} 1 & 0 \\ 0 & 1 \\ 1 & +1 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \end{array} \right). \]
Modelo, \(\mathbf{y=X\alpha+\varepsilon}\), y contrastes \(\mathbf{C^{1'}\beta}\), \(\mathbf{C^{2'}\beta}\)
Comparación entre tipo de mutantes \((\mathbf{C^{1'}\beta})\) \[\begin{array}{ccc} \beta_1^1 &=& \alpha_1-\alpha_2,\\ \beta_2^1 &=& \alpha_3-\alpha_2,\\ \beta_3^1 &=& \alpha_2-\alpha_3. \end{array}\]
Comparación entre cada mutantes y el wild type \((\mathbf{C^{2'}\beta})\) \[\begin{array}{ccc} \beta_1^2 &=& \alpha_4-\alpha_1,\\ \beta_2^2 &=& \alpha_3-\alpha_1,\\ \beta_3^2 &=& \alpha_2-\alpha_1. \end{array}\]
\[ \left( \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \\ y_7 \\ y_8 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \end{array} \right)}_{\text{Matriz de Diseño}, \mathbf{X}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right) + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6\\ \varepsilon_7\\ \varepsilon_8 \end{array} \right) \]
\[ \left( \begin{array}{r} \beta_1^1 \\ \beta_2^1 \\ \beta_3^1 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} 1 & -1 & 0 & 0\\ 1 & 0 & -1 & 0\\ 0 & 1 & -1& 0 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C^1}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right). \]
\[ \left( \begin{array}{r} \beta_1^2 \\ \beta_2^2 \\ \beta_3^2 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} 1 & 0 & 0 & -1\\ 0 & 1 & 0 & -1\\ 0 & 0 & 1 & -1 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C^2}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right). \]
El objetivo de este experimento es estudiar el efecto de las citoquinas sobre la estimulación de una sustancia (LPS) y ver como esta relación se ve afectada por la edad.
Se trata de un modelo de un un factor (tratamiento, asignable por el individuo) y un bloque (edad, no asignable, prefijada en cada ratón). En la práctica el modelo equivale al del análisis de la varianza de dos factores.
| % after : or … |
| 1 |
| 2 |
| 3 |
| 4 |
| 5 |
| 6 |
| 7 |
| 8 |
| 9 |
| 10 |
| 11 |
| 12 |
\end{center}
Aquí se adopta la segunda parametrización.
Parámetros del modelo: \[\begin{array}{ccc} \alpha_1&=& \mathbf{E} (log{LPS.Aged} ),\ \alpha_2= \mathbf{E} (log{Med.Aged} ),\\ \alpha_3&=&\mathbf{E} (log{LPS.Young} ),\, \alpha_4=\mathbf{E} (log{Med.Young} ). \end{array}\]
Contrastes: Preguntas que interesa responder: \[\begin{array}{ccc} \beta_1^1 &=& \alpha_3-\alpha_1,\quad \mbox{Efecto del tratamiento en ratones viejos} \\ \beta_2^1 &=& \alpha_4-\alpha_2,\quad \mbox{Efecto del tratamiento en ratones jóvenes} \\ \beta_3^1 &=& (\alpha_3-\alpha_1)- (\alpha_2-\alpha_4),\quad \mbox{Interacción: diferencia entre efectos} \end{array}\]
Modelo lineal:
Modelo, \(\mathbf{y=X\alpha+\varepsilon}\), y contrastes \(\mathbf{C'\beta}\)
$$ ( \[\begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \\ y_7 \\ y_8 \\ y_9 \\ y_{10} \\ y_{11} \\ y_{12} \end{array}\] ) = _{, } ( \[\begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array}\] ) + ( \[\begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6\\ \varepsilon_7\\ \varepsilon_8 \end{array}\]) $$
\[ \left( \begin{array}{r} \beta_1^1 \\ \beta_2^1 \\ \beta_3^1 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} -1 & 0 & 1 & 0\\ 0 & -1 & 0 & 1\\ -1 & 1 & 1 & -1 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right) \]
Una vez se ha expresado el experimento como un modelo lineal: \[ \mathbf{E(y_g)=X\alpha_g},\quad \text{var}(y_g)=W_g\sigma_g, \] es posible usar la teoría estándar del modelo lineal (ver (???)) para obtener:
Los procedimientos de estimación y de inferencia no dependen de qué parametrización se ha adoptado, a pesar de que distintas parametrizaciones piedan dar lugar a distintos valores numéricos..
El enfoque del modelo lineal es flexible y potente:
Por otro lado hay que tener en cuenta que, si las suposiciones no se cumplen, entonces las conclusiones deben tomarse con precaución.
Y lo que es peor, aún siendo válidas las suposiciones del modelo los resultados pueden verse afectados por el tamaño de la muestra de forma que, en muestras pequeñas donde pueden haber variaciones más grandes es fácil que se obtengan valores \(t\) no significativos o, por el contrario, excesivamente significativos (si la variación es muy reducida).
La metodología desarrollada por Smyth (Smyth, Michaud, and Scott (2005)) basada en los resultados de L"onsted & Speed (Speed (2003)) está dirigida a abordar cómo hacer frente a estas debilidades.
Smyth (Smyth, Michaud, and Scott (2005)) considera el problema de la identificación de genes que se expresan diferencialmente en las condiciones especificadas en el diseño de experimentos de microarrays. Como hemos dicho repetídamente la variabilidad de los valores de expresión difiere entre genes, pero la naturaleza paralela de la inferencia en microarrays sugiere la posibilidad de usar la información de todos los genes a la vez para mejorar la estimación de los parámetros, lo que puede llevar a una inferencia más fiable.
Básicamente lo que propone Smyth (Smyth, Michaud, and Scott (2005)) es una solución en tres pasos:
El hecho de trabajar con logaritmos permite poner el punto de corte en el cero: A mayor valor positivo más probable es que el gen esté diferencialmente expresado. A mayor valor negativo, más probable es que no lo esté
El paquete de Bioconductor limma al que hemos hecho referencia en los párrafos anteriores implementa el método de Smyth para la regularización de la varianza lo que lo ha convertido en muy popular entre los usuarios de microarrays
El código siguiente muestra como se crea la matriz del diseño y los contrastes para realizar el análisis mediante modelos lineales;
Empezamos por la matriz de diseño.
## Aged.LPS Aged.MED Young.LPS Young.MED
## Aged LPS 80L.CEL 1 0 0 0
## Aged LPS 86L.CEL 1 0 0 0
## Aged LPS 88L.CEL 1 0 0 0
## Aged Medium 81m.CEL 0 1 0 0
## Aged Medium 82m.CEL 0 1 0 0
## Aged Medium 84m.CEL 0 1 0 0
## Young LPS 75L.CEL 0 0 1 0
## Young LPS 76L.CEL 0 0 1 0
## Young LPS 77L.CEL 0 0 1 0
## Young Medium 71m.CEL 0 0 0 1
## Young Medium 72m.CEL 0 0 0 1
## Young Medium 73m.CEL 0 0 0 1
## attr(,"assign")
## [1] 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$lev
## [1] "contr.treatment"
Usando los nombres de las columnas de la matriz de diseño creamos los contrastes
require(limma)
cont.matrix <- makeContrasts (
LPS.in.AGED=(Aged.LPS-Aged.MED),
LPS.in.YOUNG=(Young.LPS-Young.MED),
AGE=(Aged.MED-Young.MED),
levels=design)
cont.matrix
## Contrasts
## Levels LPS.in.AGED LPS.in.YOUNG AGE
## Aged.LPS 1 0 0
## Aged.MED -1 0 1
## Young.LPS 0 1 0
## Young.MED 0 -1 -1
Una vez definida la matriz de diseño y los contrastes podemos pasar a estimar el modelo, estimar los contrastes y realizar las pruebas de significación que nos indiquen, para cada gen y cada comparación, si puede considerarse diferencialmente expresado.
La penúltima instrucción ejecuta el proceso de regularización de la varianza utilizando modelos de Bayes empíricos para combinar la información de toda la matriz de datos y de cada gen individual y obtener estimaciones de error mejoradas.
El análisis proporciona los estadísticos de test habituales como Fold–change \(t\)-moderados o \(p\)-valores ajustados que se utilizan para ordenar los genes de mas a menos diferencialmente expresados.
A fin de controlar el porcentaje de falsos positivos que puedan resultar del alto numero de contrastes realizados simultaneamente los p–valores se ajustan de forma que tengamos control sobre la tasa de falsos positivos utilizando el metodo de Benjamini y Hochberg (Benjamini and Hochberg (1995)).
La funcion topTable genera para cada contraste una lista de genes ordenados de mas a menos diferencialmente expresados.
topTab_LPS.in.AGED <- topTable (fit.main, number=nrow(fit.main), coef="LPS.in.AGED", adjust="fdr")
topTab_LPS.in.YOUNG <- topTable (fit.main, number=nrow(fit.main), coef="LPS.in.YOUNG", adjust="fdr")
topTab_AGE <- topTable (fit.main, number=nrow(fit.main) , coef="AGE", adjust="fdr")
Se puede visualizar los resultados mediante un volcano plot que, en este caso, representa en abscisas los cambios de expresión en escala logarítmica y en ordenadas el estadístico \(B\) en vez de el “menos logaritmo” del p-valor.
coefnum = 1
opt <- par(cex.lab = 0.7)
volcanoplot(fit.main, coef=coefnum, highlight=10, names=fit.main$ID,
main=paste("Differentially expressed genes",colnames(cont.matrix)[coefnum], sep="\n"))
abline(v=c(-1,1))
par(opt)
El análisis de datos de microarrays suele proceder, en general secuencialmente, a través de una serie de etapas tal y como se comentado en capítulos anteriores.
En cada una de las etapas pueden llevarse a cabo diversos procesos, con diversas variantes de métodos parecidos.A lo largo de la primera década del 2000 se ha desarrollado el proyecto Bioconductor que ha impulsado la creación de cientos de paquetes de R que implementan casi todas las capacidades posibles de análisis de microarrays. Aunque nadie puede conocer todos los paquetes con tan sólo conocer algunos de ellos es posible llevar a cabo potentes análisis con relativa facilidad.
El objetivo de este capítulo es ilustrar de forma breve pero completa como se hace un análisis de microarrays usando las facilidades que ofrecen Bioconductor y R.
Como se ha dicho, uno de los problemas en este tipo de estudios es que en cada etapa se puede proceder de varias formas lo que da lugar a un asfixiante número de posibilidades especialmente para el neófito. Con el fin de evitar este problema de momento se describe una sola de las opciones posibles para cada paso, lo que constituye un proceso “al estilo Bioconductor”
El ejemplo analizado consiste en unos datos disponibles en forma de paquete (estrogen) en Bioconductor. Se trata de un estudio en el que se estudia la influencia del tratamiento con estrógeno y del tiempo transcurrido desde el tratamiento en la expresión de genes asociados con cancer de mama. En capítulos anteriores se han utilizado fragmentos de este ejemplo descrito en el capítulo here en el epígrafe corresponiente a los ejemplos (ejemplo here).
El diseño experimental para este estudio consiste en un diseño factorial de 2 factores “Estrogeno” (Presente o Ausente), “Tiempo” (10h o 48h). La documentación del paquete estrogen ofrece más detalles acerca del diseño y las variables utilizadas.
Si el paquete no se encuentra instalado puede instalarse con la instrucción BiocManager::install que se descarga de la web de Bioconductor.
if (!(require("estrogen", character.only=T))){
BiocManager::install("estrogen")
}
En lo que sigue supondremos que se ha llevado a cabo una instalación estándar de Bioconductor es decir se han ejecutado las instrucciones:
BiocManager::install()
Para facilitar supondremos que trabajamos en un directorio escogido por nosotros y cuya localización se asigna a la variable workingDir. Los datos se copiaran en un subdirectorio del anterior denominado “data” que se almacenará en la variable dataDir y los resultados se almacenarán en un directorio “results” cuyo nombre completo se almacenará en la variable resultsDir.
workingDir <-getwd()
if (!file.exists("datos")) system("mkdir datos")
if (!file.exists("datos/estrogen")) system("mkdir datos/estrogen")
if (!file.exists("results")) system("mkdir results")
dataDir <-file.path(workingDir, "datos/estrogen")
resultsDir <- file.path(workingDir, "results")
setwd(workingDir)
options(width=80)
options(digits=5)
Para un análisis de datos de microarrays de Affymetrix se necesitan los archivos de imágenes escaneadas (“.CEL”) y un archivo en el que se asigne una condición experimental a cada archivo.
Una ventaja de este ejemplo es que los datos se encuentran disponibles tras instalar el paquete estrogen por lo que se pueden copiar un directorio de trabajo o simplemente trabajar con ellos en el directorio original. El directorio en que se encuentran los archivos .CEL es:
library(estrogen)
estrogenDir <- system.file("extdata", package = "estrogen")
print(estrogenDir)
## [1] "/home/alex/R/x86_64-pc-linux-gnu-library/4.0/estrogen/extdata"
Para realizar el análisis es preciso copiar todos los archivos del directorio “extdata” a nuestro directorio de trabajo “data”.
ATENCION: Esta operación que depende de la instalación y del sistema operativo. Por ejemplo para copiarlos desde R en linux se escribirá:
system(paste ("cp ", estrogenDir,"/* ./data", sep=""))
Los archivos copiados son
El proceso de leer los datos puede parecer un poco extraño a primera vista pero la idea es simple:
En primer lugar creamos algunas estructuras de datos que contienen la información sobre las variables y - opcionalmente - sobre las anotaciones y el experimento yA continuación nos basamos en estas estructuras para leer los datos y crear los objetos principales que se utilizarán para el análisis.
library(Biobase)
library(affy)
sampleInfo <- read.AnnotatedDataFrame(file.path(estrogenDir,"targLimma.txt"),
header = TRUE, row.names = 1, sep="\t")
fileNames <- pData(sampleInfo)$FileName
rawData <- read.affybatch(filenames=file.path(estrogenDir,fileNames),
phenoData=sampleInfo)
En este ejemplo se incluye un archivo defectuoso con fines didácticos, llamado “badcel”. Para ilustrar las diferencias entre archivos correctos e incorrectos se puede leer en otro objeto.
library(affy)
sampleInfoWithBad <- read.AnnotatedDataFrame(file.path(dataDir,"phenoDataWithBad.txt"),
header = TRUE, row.names = NULL, sep="\t")
fileNames <- pData(sampleInfoWithBad)$FileName
rawData.wrong <- read.affybatch(filenames=file.path(dataDir,fileNames),
phenoData=sampleInfoWithBad)
Tras leer los datos pasamos al preprocesado. Aunque puede interpretarse de distintas formas esta fase suele consistir en
La exploración previa puede hacerse paso a paso o en bloque si se utiliza algún paquete como arratQualityMetrics. En la ayuda de este paquete se encuentra una descripción de los análisis básicos que permite realizar.
El primer gráfico que suele hacerse és algun tipo de “densidad” que muestre la distribución de las señales en los datos “crudos”, es decir sin normalizar. Estos gráficos también permiten hacerse una idea de si las distribuciones de los distintos arrays son similares en forma y posición.
El gráfico de degradación –que no aparece en este caso, ya que daba problemeas al representarlo– permite hacerse una idea de como ha sido el proceso de hibridación de las muestras. Lineas paralelas sugieren una calidad similar.
## low10A low10B hi10A hi10B low48A low48B hi48A hi48B
## slope -0.103 -0.217 -0.129 -0.3810 -0.50400 -0.55300 -0.3510 -0.66700
## pvalue 0.550 0.194 0.394 0.0341 0.00482 0.00109 0.0464 0.00032
EL código para el gráfico de degradación sería:
El diagrama de cajas muestra, como el histograma, una idea de la distribución de los datos.
Finalmente un cluster jerárquico seguido de un dendrograma nos puede ayudar a hacernos una idea de si las muestras se agrupan por condiciones experimentales.
Si lo hacen es bueno, pero si no, no es necesariamente indicador de problemas, puesto que es un gráfico basado en todo los datos.
Las exploraciones anteriores nos proporcionan una idea de como son los datos.
Se pueden realizar controles de calidad más estrictos como:
El paquete arrayQualityMetrics encapsula unos cuantos análisis, de forma que con una instrucción se pueden realizar todos los análisis y enviar la salida a un directorio.
if(!require(arrayQualityMetrics)) BiocManager::install("arrayQualityMetrics")
library(arrayQualityMetrics)
arrayQualityMetrics(rawData, outdir = "Informe_de_calidad_para_los_datos_del_caso_ESTROGENO")
arrayQualityMetrics(rawDataBad, outdir = "Informe_de_calidad_(2)_para_los_datos_del_caso_ESTROGENO")
Estos informes contienen gráficos, y explicaciones de como interpretarlos. Los gráficos llevan también indicadores para determinar si una muestra determinada puede tener algún tipo de problema.
En el código de ejemplo hemos realizado dos veces el análisis, una con los datos “buenos” y otra incluyendo el array defectuoso (“bad.cel”). Como puede verse en las tablas resumen, y en los directorios de resultados el array defectuoso es detectado por la mayoría de pruebas.
La figuras ?? y ?? presenta una tabla que resumen algunas de las pruebas (no todas) realizadas e indica si, para cada criterio, puede considerarse que un array dado es un outlier.
Para concluir con esta seccion merece la pena recordar que todos estos graficos son exploratorios. Raramente se descarta un array si solo un grafico lo sugiere.
Una vez realizado el control de calidad se procede a normalizar los datos y sumarizarlos.
Hecho esto puede realizarse un filtraje no específico con el fin de eliminar genes que constituyen básicamente “ruído”, bien porque sus señales son muy bajas o bien porque apenas varían entre condiciones, por lo que no aportan nada a la selección de genes diferencialmente expresados.
La normalización puiede hacerse por distintos métodos (MAS5, VSN, RMA, GCRMA, …) En este ejemplo se utilizará el método RMA pero no se realizará filtraje alguno. Esto puede implicar quizás que para seleccionar genes diferencialmente expresados basándose en el ajuste de p-valores debamos utilizar criterios menos restrictivos que si hubiéramos filtrado, pero tiene la ventaja de eliminar un paso que, en el mejor de los casos, resulta controvertido.
El procesado mediante RMA implica un proceso en tres etapas:
library(affy)
eset_rma <- rma(rawData)
## Background correcting
## Normalizing
## Calculating Expression
Si se desea normalizar únicamente si los datos normalizados no estan disponibles puede utilizarse el código siguiente:
library(affy)
normalize <- T
if(normalize){
eset_rma <- rma(rawData)
save(eset_rma, file=file.path(dataDir,"estrogen-normalized.Rda"))
}else{
load (file=file.path(dataDir,"estrogen-normalized.Rda"))
}
La normalización hace que los valores de los arrays sean comparables entre ellos, aunque los distintos métodos situan los valores en escalas distintas, por lo que lo que no resulta directamente comparable son los valores normalizados por distintos métodos.
Un boxplot de los valores normalizados muestra que los valores ya están en una escala en donde se pueden comparar.
Ahora bien, esto no debe considerarse sinónimo de buena calidad porque si se aplica la normalización por cuantiles el gráfico aparecerá así independientemente de otros factores que puedan afectar la calidad.
boxplot(eset_rma,main="Boxplot for RMA-normalized expression values ",
names=sampleNames, cex.axis=0.7, col=info$grupo+1,las=2)
El código siguiente, que no se ejecuta muestra como se podrían aplicar distintos métodos para lego compararlos entre ellos.
eset_mas5 <- mas5(rawData) # Uses expresso (MAS 5.0 method) much slower than RMA!
stopifnot(require(gcrma))
eset_gcrma <- gcrma(rawData) # The 'library(gcrma)' needs to be loaded first.
stopifnot(require(plier))
eset_plier <- justPlier(rawData, normalize=T) # The 'library(plier)' needs to be loaded first.
compara <-data.frame(RMA=exprs(eset_rma)[,1], MAS5 =exprs(eset_mas5)[,1],
GCRMA=exprs(eset_gcrma)[,1], PLIER =exprs(eset_plier)[,1])
pairs(compara)
El filtraje no específico permite eliminar los genes que varían poco entre condiciones o que deseamos quitar por otras razones como por ejemplo que no disponemos de anotación para ellos. La función nsFilter permite eliminar los genes que, o bien varían poco, o bien no se dispone de anotación para ellos.
Si al filtrar deseamos usar las anotaciones, o la falta de ellas, como criterio de filtraje debemos disponer del correspondiente paquete de anotaciones.
require(genefilter)
filtered <- nsFilter(eset_rma, require.entrez=TRUE,
remove.dupEntrez=TRUE, var.func=IQR,
var.cutoff=0.5, var.filter=TRUE,
filterByQuantile=TRUE, feature.exclude="^AFFX")
La función nsFilter devuelve los valores filtrados en un objeto expressionSet y un informe de los resultados del filtraje.
names(filtered)
## [1] "eset" "filter.log"
class(filtered$eset)
## [1] "ExpressionSet"
## attr(,"package")
## [1] "Biobase"
print(filtered$filter.log)
## $numDupsRemoved
## [1] 2856
##
## $numLowVar
## [1] 4292
##
## $numRemoved.ENTREZID
## [1] 1166
##
## $feature.exclude
## [1] 19
eset_filtered <-filtered$eset
Podemos grabar el objeto eset_rma y los datos filtrados para su posterior uso.
save(eset_rma, eset_filtered, file=file.path(resultsDir, "estrogen-normalized.Rda"))
Despues del filtraje han quedado 4409 genes disponibles para analizar.
Como en las etapas anteriores la selección de genes diferencialmente expresados (GDE) puede basarse en distintas aproximaciones, desde la \(t\) de Student al programa SAM pasando por multitud de variantes.
En este ejemplo se aplicará la aproximación presentada por Smyth et al. (2004) basado en la utilización del combinada con un método para obtener una estimación mejorada de la varianza.
El capítulo 6 de tse manual o el manual del programa limma contienen explicaciones detalladas sobre construir un modelo lineal para este problema y como utilizarlo para seleccionar genes diferencialmente expresados.
El primer paso para el análisis es crear la matriz de diseño.
La situación discutida en este ejemplo se puede modelizar de dos formas, tal como se discute en la presentación citada más arriba:
Tal como se describe en el manual de limma la segunda parametrización resulta a menudo más comoda, a pesar de parecer menos intuitiva, porque permite formular con más facilidad que la de dos factores las preguntas que típicamente interesan a los investigadores.
El modelo lineal para este estudio será:
\[ \left( \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \\ y_7 \\ y_8 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \end{array} \right)}_{\mbox{Design Matrix}, \mathbf{X}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right) + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6\\ \varepsilon_7\\ \varepsilon_8 \end{array} \right) \]
Los parámetros del modelos representan las cuatro combinaciones tiempo/estrogeno.
\[\begin{array}{ccc} \alpha_1&=& \mathbf{E} (log{Abs.10h} ),\\ \alpha_2&=& \mathbf{E} (log{Pres.10h} ),\\ \alpha_3&=& \mathbf{E} (log{Abs.48h} ),\\ \alpha_4&=& \mathbf{E} (log{Pres.48h} ). \end{array}\]
La matriz de diseño puede definirse manualmente o a partir de un factor creado específicamente para ello.
Manualmente, seria:
design.1<-matrix(
c(1,1,0,0,0,0,0,0,
0,0,1,1,0,0,0,0,
0,0,0,0,1,1,0,0,
0,0,0,0,0,0,1,1),
nrow=8,
byrow=F)
colnames(design.1)<-c("neg10h", "est10h", "neg48h", "est48h")
rownames(design.1) <- c("low10A", "low10B", "hi10A" , "hi10B", "low48A", "low48B", "hi48A" , "hi48B")
print(design.1)
## neg10h est10h neg48h est48h
## low10A 1 0 0 0
## low10B 1 0 0 0
## hi10A 0 1 0 0
## hi10B 0 1 0 0
## low48A 0 0 1 0
## low48B 0 0 1 0
## hi48A 0 0 0 1
## hi48B 0 0 0 1
Alternativamente puede crearse la matriz de diseño a partir de la información sobre las condiciones contenida en el phenoData, siempre que exista un campo adecuado para ello.
En este caso la columna Target se ha creado para utilizarla con esta finalidad. El objeto phenoData puede recrearse a partir del archivo original o extrayéndolo del objeto ExpresionSet que contiene los datos y las covariables.
## neg10h est10h neg48h est48h
## low10A 1 0 0 0
## low10B 1 0 0 0
## hi10A 0 1 0 0
## hi10B 0 1 0 0
## low48A 0 0 1 0
## low48B 0 0 1 0
## hi48A 0 0 0 1
## hi48B 0 0 0 1
## attr(,"assign")
## [1] 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$lev
## [1] "contr.treatment"
Ambas matrices, design, design.1 resultan iguales.
Dado un modelo lineal definido a traves de una matriz de diseño pueden formularse las preguntas de interes como contrastes es decir comparaciones entre los parametros del modelo.
Cada parametrizacion distinta requerira de unos contrastes diferentes para las mismas preguntas, por lo que habitualmente se utilizara la parametrizacion que permita formular de forma mas clara las comparaciones de interes.
En este caso interesa estudiar
Esto se puede formular facilmente con la parametrizacion adoptada. \[\begin{array}{ccc} \beta_1^1 &=& \alpha_2-\alpha_1,\quad \mbox{Efecto del estrogeno pasadas 10 horas} \\ \beta_2^1 &=& \alpha_4-\alpha_3,\quad \mbox{Efecto del estrogeno pasadas 48 horas} \\ \beta_3^1 &=& \alpha_3-\alpha_1,\quad \mbox{Efecto del tiempo en ausencia de Estrogeno} \end{array}\]
cont.matrix <- makeContrasts (
Estro10=(est10h-neg10h),
Estro48=(est48h-neg48h),
Tiempo=(neg48h-neg10h),
levels=design)
cont.matrix
## Contrasts
## Levels Estro10 Estro48 Tiempo
## neg10h -1 0 -1
## est10h 1 0 0
## neg48h 0 -1 1
## est48h 0 1 0
Una vez definida la matriz de diseño y los contrastes podemos pasar a estimar el modelo, estimar los contrastes y realizar las pruebas de significación que nos indiquen, para cada gen y cada comparación, si puede considerarse diferencialmente expresado.
require(limma)
fit<-lmFit(eset_rma, design)
fit.main<-contrasts.fit(fit, cont.matrix)
fit.main<-eBayes(fit.main)
El método implementado en amplía el análisis tradicional utilizando modelos de Bayes empíricos para combinar la información de toda la matriz de datos y de cada gen individual y obtener estimaciones de error mejoradas.
El análisis proporciona los estadísticos de test habituales como Fold–change \(t\)-moderados o \(p\)-valores ajustados que se utilizan para ordenar los genes de mas a menos diferencialmente expresados.
A fin de controlar el porcentaje de falsos positivos que puedan resultar del alto numero de contrastes realizados simultaneamente los p–valores se ajustan de forma que tengamos control sobre la tasa de falsos positivos utilizando el metodo de Benjamini y Hochberg (Benjamini and Hochberg (1995)).
La funcion topTable genera para cada contraste una lista de genes ordenados de mas a menos diferencialmente expresados.
topTabEstro10 <- topTable (fit.main, number=nrow(fit.main), coef="Estro10", adjust="fdr")
topTabEstro48 <- topTable (fit.main, number=nrow(fit.main), coef="Estro48", adjust="fdr")
topTabTiempo <- topTable (fit.main, number=nrow(fit.main) , coef="Tiempo", adjust="fdr")
Una forma de visualizar los resultados es mediante un volcano plot que representa en abscisas los cambios de expresión en escala logarítmica y en ordenadas el “menos logaritmo” del p-valor o alternativamente el estadístico \(B\).
Cuando se realizan varias comparaciones a la vez puede resultar importante ver que genes cambian simultáneamente en más de una comparación. Si el número de comparaciones es alto también puede ser necesario realizar un ajuste de p-valores entre las comparaciones, distinto del realizado entre genes.
La función decidetests permite realizar ambas cosas. En este ejemplo no se ajustaran los p–valores entre comparaciones. Tan solo se seleccionaran los genes que cambian en una o más condiciones.
EL resultado del análisis es una tabla res que para cada gen y cada comparación contiene un 1 (si el gen esta sobre-expresado o “up” en esta condicion), un 0 (si no hay cambio significativo) o un -1 (si esta “down”-regulado).
Para resumir dicho análisis podemos contar qué filas tienen como mínimo una celda distinta de cero:
sum.res.rows<-apply(abs(res),1,sum)
res.selected<-res[sum.res.rows!=0,]
print(summary(res))
## Estro10 Estro48 Tiempo
## Down 23 99 27
## NotSig 12512 12375 12591
## Up 90 151 7
Un diagrama de Venn permite visualizar la tabla anterior sin diferenciar entre genes “up” o “down” regulados.
\begin{figure}
vennDiagram (res.selected[,1:3], main="Genes in common #1", cex=0.9)
\end{figure}
Alizadeh, A., M. B. Eisen, E. Davis, C. Ma, I. Lossos, A. Rosenwald, J. Boldrick, et al. 2000. “Distinct Types of Diffuse Large B–Cell Lymphoma Identified by Gene Expression Profiling.” Nature 403: 503–11.
Allison, David B, Xiangqin Cui, Grier P Page, and Mahyar Sabripour. 2006. “Microarray Data Analysis: From Disarray to Consolidation and Consensus.” Nat Rev Genet 7 (1): 55–65.
Benjamini, Y., and Y. Hochberg. 1995. “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” Journal of the Royal Statistical Society B 57: 289–300.
Chelvarajan, R. L., Y. Liu, D. Popa, M. L. Getchell, T. V. Getchell, A. J. Stromberg, and S. Bondada. 2006. “Molecular basis of age-associated cytokine dysregulation in LPS-stimulated macrophages.” Journal of Leukocyte Biology 79 (6): 1314.
Dudoit, S., J. P. Shaffer, and J. C. Boldrick. 2003. “Multiple Hypothesis Testing in Microarray Experiments.” Statistical Science 18: 71–103.
Golub, T. R., D. K. Slonim, P. Tamayo, C. Huard, M. Gaasenbeek, J. P. Mesirov, H. Coller, et al. 1999. “Molecular Classification of Cancer: Class Discovery and Class Prediction by Gene Expression Monitoring.” Science 286: 531–37.
Imbeaud, Sandrine, and Charles Auffray. 2005. “’The 39 Steps’ in Gene Expression Profiling: Critical Issues and Proposed Best Practices for Microarray Experiments.” Drug Discovery Today 10 (17): 1175–82. https://doi.org/10.1016/S1359-6446(05)03565-8.
Kerr, M K, and G A Churchill. 2001. “Experimental Design for Gene Expression Microarrays.” Biostatistics.
Moore, Shanna, Julia Vrebalov, Paxton Payton, and Jim Giovannoni. 2002. “Use of Genomics Tools to Isolate Key Ripening Genes and Analyse Fruit Maturation in Tomato.” Journal of Experimental Botany 53 (377): 2023–30. https://doi.org/10.1093/jxb/erf057.
Quackenbush, J. 2002. “Microarray Data Normalization and Transformation.” Nature Genet. 32: 496–501.
Scholtens, Denise, Alexander Miron, Faisal M. Merchant, Arden Miller, Penelope L. Miron, J. Dirk Iglehart, and Robert Gentleman. 2004. “Analyzing Factorial Designed Microarray Experiments.” J. Multivar. Anal. 90 (1): 19–43. https://doi.org/http://dx.doi.org/10.1016/j.jmva.2004.02.004.
Simon, Richard M., Edward L. Korn, Lisa M. McShane, Michael D. Radmacher, George W. Wright, and Yingdong Zhao. 2003. Design and Analysis of Dna Microarray Investigations. Springer-Verlag.
Smyth, Gordon K, Joëlle Michaud, and Hamish S Scott. 2005. “Use of Within-Array Replicate Spots for Assessing Differential Expression in Microarray Experiments.” Bioinformatics 21 (9): 2067–75.
Speed, T. 2003. Statistical Analysis of Gene Expression Data. Boca Raton, Fla.: Chapman & Hall/CRC.
Tibshirani, R. 2006. “A simple method for assessing sample sizes in microarray experiments.” BMC Bioinformatics 7 (1): 106.
van ’t Veer, Laura J, Hongyue Dai, Marc J van de Vijver, Yudong D He, Augustinus A M Hart, Mao Mao, Hans L Peterse, et al. 2002. “Gene Expression Profiling Predicts Clinical Outcome of Breast Cancer.” Nature 415 (6871): 530–6.
Yang, Y. H., M. J. Buckley, S. Dudoit, and T. P. Speed. 2002. “Comparison of Methods for Image Analysis on cDNA Microarray Data.” Journal of Computational and Graphical Statistic 11 (1).
Zhu, Qianqian, Jeffrey Miecznikowski, and Marc Halfon. 2010. “Preferred Analysis Methods for Affymetrix Genechips. II. An Expanded, Balanced, Wholly-Defined Spike-in Dataset.” BMC Bioinformatics 11 (1): 285. https://doi.org/10.1186/1471-2105-11-285.